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Abstract:  Signals  from  many  military  targets  of  interest  are  often 
strongly  randomized,  due  to  the  irregular  mechanisms  by  which  the 
signals  are  generated  and  propagated.  In  particular,  complicated  and 
dynamic  terrestrial/atmospheric  environments  (with  man-made  objects, 
vegetation,  and  turbulence)  randomize  signals  through  random 
atmospheric  and  terrestrial  processes  affecting  the  propagation.  Signals 
may  also  be  considered  random  due  to  uncertainties  in  the  knowledge  of 
the  propagation  environment  and  target-sensor  geometry.  Predictions  of 
sensor  performance  and  recommendations  of  sensor  types  and  placements 
derived  from  them,  thus,  should  account  for  the  random  nature  of  the 
sensed  signals.  This  report  discusses  software-modeling  approaches  for 
characterizing  signals  subject  to  random  generation  and  propagation 
mechanisms.  By  representing  signals  with  random  variables,  they  are 
manipulated  statistically  to  make  probabilistic  predictions  of  sensor 
performance.  Both  the  theory  and  implementation  in  a  general,  object- 
oriented  software  design  for  battlefield  signal  transmission  and  sensing  is 
explained.  The  Java-language  software  program  is  called  Environmental 
Awareness  for  Sensor  and  Emitter  Employment.  Some  important 
numerical  issues  in  the  implementation  are  also  discussed. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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1  Introduction 

A  variety  of  modern-day  Army  missions  rely  on  effective  sensing  capa¬ 
bilities  that  provide  intelligence  on  the  adversary  and  to  protect  friendly 
forces  from  enemy  detection.  Sensors  that  are  stationary  (e.g.,  micro¬ 
phones,  geophones,  and  ground-based  radars)  and  moving  (e.g.,  cameras 
on  unattended  aerial  vehicles  and  ground  vehicles)  assist  operations  such 
as  persistent  surveillance  of  small,  forward-operating  bases  and  rapid 
covert  troop  maneuvers  in  the  air  and  on  the  ground.  When  advantageous, 
sensing  is  often  performed  in  multiple  signal  modalities  including  visible, 
infrared,  acoustic,  seismic,  radiofrequency,  chemical,  and  biological. 

Yet,  despite  the  increasing  assortment  of  sensors  available,  knowledge  and 
expertise  of  how  to  use  them  efficiently  in  particular  environments  and 
missions  is  quite  frequently  lacking  (Hieb  et  al.  2007).  Since  terrain  and 
weather  effects  on  signals  are  complex  and  oftentimes  counterintuitive, 
computational  simulations  are  valuable  to  both  facilitate  quick  and  accu¬ 
rate  decision  making  when  planning  the  use  of  sensors  in  an  actual  mis¬ 
sion  as  well  as  to  guide  the  development  of  sensor  technologies  in  general. 
When  multiple  diverse  signal  modalities  are  involved,  it  would  be  ideal  to 
integrate  them  all  into  a  single  simulation  environment  rather  than  having 
one  for  each  modality. 

Environmental  Awareness  of  Sensor  and  Emitter  Employment  (EASEE)  is 
a  software  framework  that  provides  a  single  environment  for  analyzing 
sensor  performance  involving  many  different  signal  modalities.  This  abili¬ 
ty  for  multimodal  signal  analysis  then  enables  EASEE  to  also  perform 
higher-level  data  synthesis  needed  to  answer  critical  sensing  questions 
such  as  the  sensor  types  and  locations  best  suited  for  accomplishing  mis¬ 
sion  objectives  within  the  constraints  of  a  particular  environment.  A  com¬ 
prehensive  description  of  the  software  design  and  structure  of  EASEE  is 
provided  in  Wilson  et  al.  (2009). 

In  addition  to  reconciling  models  of  many,  differing  signal  modalities  into 
a  common  software  framework,  a  decision-support  tool  (DST)  must  also 
systematically  assess  the  uncertainty  with  its  predictions.  Accounting  for 
and  measuring  uncertainty  is  often  overlooked  in  current  software  tools 
predicting  sensor  performance,  but  the  reliability  of  a  DST’s  recommenda- 
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tions  is  obviously  key  to  well-informed  decision  making.  In  the  compli¬ 
cated  process  of  modeling  environmental  effects  on  signal  transmission 
and  sensing,  many  types  and  levels  of  uncertainty  are  involved.  A  thorough 
explanation  of  practical  uncertainty  issues  is  available  in  Wilson  et  al. 
(2008). 

Generally,  uncertainties  may  be  due  to  incomplete  and/or  inaccurate 
knowledge  (of,  say,  the  weather  for  a  particular  time  and  location)  or  due 
to  processes  that  are  purely  stochastic  (e.g.,  sensitivity  of  wave  propaga¬ 
tion  to  irresolvable  and  highly  variable  environmental  details).  It  is  the  lat¬ 
ter  type  of  uncertainty— namely  stochastic/random  uncertainty— with 
which  this  report  is  concerned.  (Note  that  sensor  performance  is  also  af¬ 
fected  by  non-stochastic  processes  not  addressed  here.)  In  principle,  prop¬ 
agation  of  acoustic,  seismic,  and  electromagnetic  (including  visual  and 
infrared)  waves  as  well  as  the  dispersion  of  transient  gases  are  still  sensi¬ 
tive  to  many  fine-scale  and  even  dynamic  environmental  variations  that 
cannot  be  known  or  determined.  Practically  speaking,  second-to-second 
variability  in  signal  characteristics  (e.g.,  sound  level,  seismic  energy,  con¬ 
centration  of  a  chemical  agent,  etc.)  caused  by  unobservable  and  unpre¬ 
dictable  environmental  features  (e.g.,  turbulent  eddies,  vegetation,  dust 
particles,  precipitation,  small  urban  structures)  are  essentially  random,  so 
statistical  models  are  necessary  to  represent  signal  and  noise  (Strohbehn 
1978;  Andrews  and  Phillips  2005).  Wilson  et  al.  (2008)  outlines  specific 
examples  in  greater  detail.  Probabilistic  information  about  signal  and 
noise  may  then  be  processed  to  compute  statistical  metrics  on  sensor  per¬ 
formance  and  signal-emitter  detectability. 

The  report  begins  by  first  explaining  the  premise  and  theory  of  modeling 
signals  statistically  in  Chapter  2,  which  also  describes  how  statistical  mod¬ 
eling  of  signals  lead  to  probabilistic  predictions  on  sensor  performance. 
Then,  Chapter  3  explains  how  the  statistical  modeling  is  implemented  effi¬ 
ciently  in  EASEE’s  general,  object-oriented  structure.  A  discussion  on  var¬ 
ious  approaches  for  summing  multiple  random  variables  is  given  in  Chap¬ 
ter  4,  which  is  required  when  simulating  multiple  target  and  noise  signals 
arriving  at  a  sensor.  Some  general  numerical  issues  in  programming  statis¬ 
tical  computations  are  illustrated  in  Chapter  5,  while  the  Appendix  sup¬ 
plies  specific  information  on  how  particular  statistical  calculations  were 
implemented. 
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2  Statistical  Modeling  of  Signals 

Signature  data  features 

The  EASEE  software  design  (Wilson  et  al.  2009)  operates  on  basic,  identi¬ 
fying  qualities  of  a  signal  called  signature  data  features  (or  features  for 
short).  The  features  can  be  thought  of  as  the  most  essential  characteristics 
of  signals  that  might  identify  their  source  (i.e.,  an  observable  that  is  a  func¬ 
tion  of  time  at  the  source,  sensor,  or  somewhere  in  between).  Generally, 
features  are  extracted  from  raw  sensor  data  after  some  low-level 
processing  (e.g.,  calibrations  to  remove  sensor  response  or  filtering  into 
spectral  bands),  typically  being  a  conservative  quantity  (such  as  energy  or 
mass).  Examples  include:  sound  power  of  a  harmonic  line  or  in  a  standard 
octave  band;  infrared  brightness  in  near,  shortwave,  midwave,  longwave, 
or  far  bands;  components  of  an  electromagnetic  field  vector;  and  concen¬ 
tration  of  a  particular  chemical  or  biological  species.  By  using  features  ra¬ 
ther  than  raw  signals,  simulations  like  EASEE  are  made  not  only  quicker 
and  much  more  efficient  but  also  more  general  and  nonspecific  to  any  par¬ 
ticular  sensor  system. 

Probabilistic  modeling  of  features 

As  discussed  earlier,  irregular  generation  mechanisms  and  random  propa¬ 
gation  processes  may  make  it  impractical  to  predict  feature  characteristics 
deterministically;  rather,  they  can  only  be  described  probabilistically . 
Thus,  signal  and  noise  features  may  be  regarded  as  random  variables. 
Then,  statistical  distributions  of  signal  features  are  generated,  propagated, 
and  processed. 

A  variety  of  statistical  models  are  appropriate  for  representing  sensor  data. 
For  example,  the  exponential  distribution  describes  a  single,  strongly  scat¬ 
tered  signal,  which  often  arises  when  acoustic  or  electromagnetic  waves 
are  scattered  by  both  objects  and  turbulent  wind.  A  version  of  the  Rice- 
Nakagami  model  (specifically  with  a  variable  transformation  of  x  X2) 
may  closely  approximate  signal  power  distributions  of  certain  acoustic  and 
seismic  signals.  Since  the  lognormal  distribution  models  variables  that  can 
be  the  multiplicative  product  of  many  independent,  positive  random  va¬ 
riables,  atmospheric  scientists  often  use  it  to  describe  plume  sizes  and  fre¬ 
quency  distributions  of  transient  gases  from  turbulent  processes  in  the  air 
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(Baker  et  al.  1983;  Limpert  et  al.  2001).  Generally,  signals  are  represented 
with  probability  density  functions  (pdfs)  that  have  nonnegative  support 
because  signal  power  or  concentration  can  never  be  negative.  However,  if 
the  mean  is  positive  and  much  larger  than  the  standard  deviation,  pdfs 
with  real-number  support  like  the  Gaussian  may  also  be  suitable  for 
representing  a  signal.  Currently,  EASEE  implements  five  continuous  sta¬ 
tistical  models— namely,  the  Gaussian,  lognormal,  exponential,  gamma, 
and  the  jr  x2  transformed  Rice-Nakagami— as  well  as  a  discrete  model. 
(Other  examples  of  statistical  models  for  signals  are  in  Burdic  1984;  Wil¬ 
son  et.  al.  2002;  and  Nadarajah  2008) 

For  consistency,  all  of  these  models  are  used  to  describe  distributions  of 
signal  power  (or  intensity)  rather  than  amplitude.  Thus,  while  the  untrans¬ 
formed  Rice-Nakagami  model  may  describe  the  distribution  of  signal  am¬ 
plitudes  of  various  acoustic  and  seismic  signals,  it  is  necessary  to  find  the 
corresponding  model  that  would  represent  the  distribution  of  the  signal 
powers.  Since  signal  power  is  related  to  the  square  of  the  amplitude,  the 
correct  model  for  the  signal-power  distribution  is  the  Rice-Nakagami  with 
a  variable  transformation  of  X  X2 ,  which  is  what  is  coded  in  EASEE. 
(Whenever  the  “transformed  Rice-Nakagami”  is  referred  to  in  this  report, 
the  X->X2  transformed  Rice-Nakagami  model  is  meant.)  Details,  includ¬ 
ing  the  derivation  of  the  X->X2  transformed  Rice-Nakagami  random  vari¬ 
able,  are  provided  in  the  Appendix. 

Statistical  analysis  for  measuring  sensor  performance 

When  features  are  described  probabilistically,  it  is  subsequently  possible 
to  derive  statistical  metrics  of  sensor  performance.  There  are  actually  four 
probabilities  that  are  generally  of  interest: 

1.  Fed  =  probability  that  noise  only  is  present  when  the  sensor  determines 
that  noise  only  is  present  (a  correct  dismissal) 

2.  Pfa  =  probability  that  noise  only  is  present  when  the  sensor  determines  that 
a  target  is  present  (a  false  alarm) 

3.  Pfd  =  probability  that  a  target  is  present  when  the  sensor  determines  that 
noise  only  is  present  (a  false  dismissal) 

4.  Pd  =  probability  that  a  target  is  present  when  the  sensor  determines  that  a 
target  is  present  (a  correct  detection) 
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These  probabilities  are  given  by  the  following  integrals: 

=  \f0{x)dx  =  \-P(a 
0 

oo 

pa  =lfo(x)dx  =  l-Pc, 

P 

P 

pm  =  j  f  (x)dx  =  I  -  P(i 
0 

co 

PA  =\fi{x)dx  =  l-Pa 
P 

where:  f0  =  pdf  of  the  noise  signal  alone, 

fi  =  pdf  of  the  target  and  noise  signals  together,  and 

ft  =  the  detection  threshold  set  by  the  detection  algorithm. 

A  detection  algorithm  computes  an  appropriate  [i  that  ideally  optimizes 
sensor  performance  by  increasing  probability  of  detection  and  decreasing 
probability  of  false  alarm.  Some  examples  (which  are  currently  imple¬ 
mented  in  EASEE)  include: 

1.  Neyman -Pearson  (constant  false-alarm  rate)  criterion 

2.  absolute  threshold  detection 

3.  relative  threshold  detection 

4.  error  minimization 

5.  Bayes  risk  minimization 

As  the  integrals  above  suggest,  various  probabilities  of  interest  are  com¬ 
puted  by  manipulating  pdfs  of  signal  and  noise.  In  fact,  integrations  of 
pdfs  below  and  above  a  threshold  value  are  the  cumulative  distribution 
function  (cdf)  and  the  complementary  cumulative  distribution  function 
(ccdf),  respectively.  Thus,  we  may  also  write: 

Pca  =  \f0{x)dx  =  F0(j3) 

0 

00 
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where:  /  =  pdf, 

F  =  cdf,  and 
F c  =  ccdf 

with  the  subscripts  O  and  l  denoting  the  pdf,  cdf,  or  ccdf  of  the  noise  signal 
and  combined  target  and  noise  signal,  respectively.  Strictly  speaking,  the 
cdf  is  actually  the  integral  of  the  pdf  from  negative  infinity  to  the  threshold 
value;  however,  the  integral  from  zero  would  only  be  negligibly  different  (if 
at  all)  when  signals  are  represented  with  appropriate  statistical  models, 
whose  measure  limits  to  zero  (or  is  even  nonexistent)  in  the  negative  do¬ 
main,  since  negative  signal  power  or  concentration  is  impossible.  The 
aforementioned  detection  algorithms  also  compute  the  quantile  (or  in- 
verse-cdf)  function  when  determining  signal-power  (or  concentration) 
thresholds  corresponding  to  specified  probabilities. 

To  obtain  a  particular  statistical  representation,  both  the  type  of  random 
variable  and  the  values  for  its  parameters  are  needed.  While  the  pdf  and 
cdf  for  each  type  of  random  variable  are  coded  to  use  its  unique  parame¬ 
ters  (e.g.,  location,  shape,  scale,  etc.),  it  is  useful  to  also  describe  random 
variables  with  the  universal  parameters,  mean  and  variance.  For  two- 
parameter  random  variables,  closed-form  expressions  may  be  found  to 
convert  its  unique  parameters  to/from  mean  and  variance. 

Finally,  the  signal  pdfs  within  the  integrals  above  are  often  joint  pdfs,  since 
multiple  signals  of  interest  and  noise  may  arrive  at  a  sensor  simultaneous¬ 
ly.  When  multiple  signal  sources  are  present  within  a  sensor’s  vicinity, 
random  variables  representing  each  signal  source  must  be  summed  at  the 
receiver.  Depending  on  the  type  of  random  variables  being  summed,  com¬ 
puting  the  exact  summation  can  be  complicated  and  intensive.  Thus,  ob¬ 
taining  a  quick  approximation  may  be  more  practical.  A  couple  of  approx¬ 
imating  approaches  and  the  most  efficient  method  for  computing  exact 
results  are  described  later  in  this  report. 
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3  Software  Implementation 

The  EASEE  software  design  is  formulated  within  the  conceptual  frame¬ 
work  of  object-oriented  programming  in  the  Java  language.  For  a  compre¬ 
hensive  description  on  the  object-oriented  structure  and  approach  to  sig¬ 
nal  transmission  and  sensing  in  EASEE,  it  is  recommended  to  consult 
Wilson  et  al.  (2009).  Here,  only  a  more  thorough  and  updated  explanation 
of  the  signal  model  objects  will  be  given. 

Signal  model  objects 

As  previously  described,  random  environmental  effects  on  transmitted  and 
received  features  are  accounted  for  by  representing  them  as  random 
variables.  Parametric  descriptions  of  these  random  variables  are  pro¬ 
grammed  in  EASEE  within  Java  classes.  These  Java  classes  are  denoted  as 
signal  models  and  instances  of  them  are  signal  model  objects.  Signal 
model  objects  are  key  in  the  architecture  of  EASEE,  since  they  are  what  are 
actually  transmitted,  received,  and  processed,  as  outlined  in  Figure  1.  The 
purple  arrows  in  the  figure  specify  where  signal  models  are  passed  from 
one  of  the  five  components  (illustrated  in  boxes)  to  another.  Specifically, 
th e  feature  generator  produces  a  signal  model  object  that  is  inputted  and 
then  altered  by  th  e  feature  propagator  and  feature  sensor  before  it  is 
finally  analyzed  by  th e  feature  processor  to  produce  an  inference,  which 
represents  desired  information  derived  from  the  data  features  such  as 
probability  of  detection  or  error  of  target  location  estimates. 

Signal  model  classes  contain  methods  for  the  various,  previously  described 
statistical  operations  necessary  for  making  probabilistic  predictions  on 
sensor  performance,  including: 

1.  setting  the  mean  and  variance 

2.  converting  from  mean  and  variance  to  unique  parameter  values 

3.  computing  the  pdf,  cdf,  and  quantile 

4.  summing  a  random  variable  described  by  an  instance  of  the  class  with 
another  one. 
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Figure  1.  Generalized  flow  diagram  for  transmitting,  receiving,  and  processing  signature  data 

features  in  EASEE. 


These  signal  model  objects  also  store  arrays  of  parameter  values  defining 
different  random  variables  of  the  same  type,  which  may  later  be  distri¬ 
buted  across  a  terrain  grid.  Computations  for  each  array  component  are 
parallelized  when  possible  to  improve  computational  efficiency.  As  men¬ 
tioned  earlier,  EASEE  currently  has  implementations  for  five  continuous 
statistical  models— namely,  the  Gaussian,  lognormal,  exponential,  gamma 
and  the  x^x2  transformed  Rice-Nakagami— as  well  as  a  discrete  model. 
Each  of  these  models  is  represented  via  its  own  signal  model  class. 

Signal  model  inheritance  relationships 

The  overall  signal-transmission  and  sensing  process  in  EASEE  is  described 
generally  enough  to  apply  for  all  signal  modalities  described  by  potentially 
very  different  statistics  so  that  each  of  the  generalized  components  in  the 
above  figure  is  coded  to  work  with  an  abstract  representation  of  all  statis¬ 
tical  models.  This  abstract  representation  is  the  parent  abstract  Java  class 
of  all  signal  models  and  is  called  the  abstract  signal  model.  It  defines  ab¬ 
stract  methods  for  the  various  required  statistical  operations,  which  are 
then  individually  implemented  or  overridden  by  subclass  signal  models 
representing  different  statistical  models.  Figure  2  shows  the  inheritance 
relationships  of  the  signal  models  currently  available  in  EASEE. 


Figure  2.  Inheritance  tree  of  statistical  signal  models. 
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The  constant  signal  model  that  extends  the  base  abstract  signal  model 
represent  constant  (i.e.,  time-invariant)  signals  involving  only  one  variable 
for  the  signal  mean.  A  collection  of  variable  signal  models  extend  the  con¬ 
stant  signal.  Presently,  the  inheritance  relationships  are  organized  so  that 
deeper  subclasses  add  more  characteristics  to  their  parent  classes.  Sub¬ 
classes  “inherit”  and  then  modify  the  fields  and/or  methods  of  its  parent 
class.  For  example,  nonconstant  signal  models  allow  nonzero  variance; 
whereas,  the  constant  signal  model  does  not.  The  gamma  distribution  ge¬ 
neralizes  the  exponential  distribution,  which  is,  in  fact,  a  special  case  of 
the  gamma  distribution  with  shape  parameter  equal  to  l  and  the  scale  pa¬ 
rameter  equal  to  the  mean. 

It  may  be  useful,  however,  to  structure  inheritance  relationships  different¬ 
ly.  For  example,  it  may  be  useful  to  somehow  categorize  various  distribu¬ 
tions  into  similar  statistical  families  so  the  sum  of  multiple  closely  related 
random  variables  may  be  approximated  by  a  parametric  pdf  within  the 
same  family.  One  of  the  simplest  methods  for  summing  multiple  random 
variables  (described  later)  may,  in  fact,  be  made  more  efficient  if  inherit¬ 
ance  is  organized  in  such  a  way.  Otherwise,  it  may  be  more  appropriate  to 
use  interfaces  to  delineate  statistical  families  within  the  current  signal 
model  inheritance  tree,  which  can  be  used  to  make  the  described  approxi¬ 
mation  method  for  summing  random  variables  more  efficient. 

While  the  generalized  components  in  the  simulation  scheme  are  coded  to 
accept  and  return  abstract  signal  models,  specific  subclass  signal  models 
will  actually  be  used  in  actual  simulations  as  appropriate  via  polymor¬ 
phism.  This  programming  approach  using  inheritance  and  polymorphism 
greatly  improves  software  efficiency  and  the  potential  for  EASEE  to  conti¬ 
nuously  develop  and  include  new  capabilities.  Namely,  existing  signal- 
generation,  propagation,  and  processing  algorithms  would  automatically 
operate  on  all  new  signal  models  that  extend  the  abstract  signal  model 
and,  conversely,  new  processing  methods  that  work  with  the  abstract  sig¬ 
nal  model  would  instantly  apply  to  all  existing  and  future  subclass  signal 
models. 

Representing  a  specific  feature  via  a  signal  model 

While  the  signal  models  fully  describe  a  parametric  statistical  model,  they 
are  not  explicitly  associated  with  any  particular  signature  data  feature. 
Features  are  separately  defined  as  Java -enumerated  types,  whose  repre¬ 
sentation  by  a  signal  model  are  assigned  by  a  given  subclass  implementa- 
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tion  of  the  feature  generator  as  appropriate.  Thus,  the  “signal  data  trans¬ 
missions”  passed  from  one  generalized  component  to  another  in  Figure  l 
are  essentially  “packets”  of  information  assembled  by  the  feature  genera¬ 
tor  that  lists  transmitted  features  and  their  corresponding  signal-model 
representations.  For  organization  and  utility,  features  are  defined  and 
grouped  in  different  Java  enum  classes,  typically  by  modality.  For  in¬ 
stance,  all  seismic  features  and  all  radiofrequency  features  are  enumerated 
in  their  own,  separate  classes.  When  appropriate,  further  distinctions  are 
defined  within  a  single  signal  modality  such  as  with  acoustic  features, 
which  have  separate  enum  classes  for  acoustic  octave  bands  (both  stan¬ 
dard  and  third-octave)  and  linearly  spaced  spectral  bands. 


ERDC/CRREL  TR-10-12 


11 


4  Sums  of  Random  Variables 

When  simulating  multiple  signals  of  interest  and  noise  arriving  at  a  sen¬ 
sor,  probability-of-detection  calculations  are  performed  on  the  pdf  of  the 
sum  of  multiple  signals.  Thus,  the  final  pdf  that  is  analyzed  at  the  sensor 
represents  a  sum  of  multiple  random  variables,  where  each  individual  sig¬ 
nal  is  described  by  a  unique  pdf.  The  calculation  is  trivial  for  summing 
multiple  Gaussian  random  variables,  since  the  pdf  of  the  sum  is  another 
Gaussian  pdf  with  mean  and  variance  equal  to  the  sum  of  the  means  and 
variances,  respectively,  of  each  contributing  random  variable.  In  fact,  the 
mean  and  variance  of  the  sum  of  any  (even  potentially  dissimilar)  types  of 
random  variables  is  always,  by  definition,  the  sum  of  the  means  and  va¬ 
riances  of  each  summed  random  variable.  (This  property  is  referred  to  in 
the  rest  of  this  report  by  saying  that  the  means  and  variances  of  the 
summed  random  variables  are  “conserved.”)  Yet,  despite  this  property, 
implementation  for  summing  multiple  non-Gaussian  random  variables  is 
not  straightforward  because  an  analytical  form  for  the  pdf  of  the  sum  is 
usually  not  available. 

A  very  simple  scheme 

One  of  the  simplest  methods  for  representing  the  sum  of  multiple  random 
variables  of  related  types  is  to  approximate  the  sum  with  a  pdf  that  con¬ 
serves  the  mean  and  variance  of  the  original  contributions.  The  pdf  for  the 
sum  may  be  the  same  as  one  of  the  original  random  variables,  or  a  genera¬ 
lized  form  of  the  pdf  of  the  original  variables.  This  approach  is  exact  if  all 
the  summed  random  variables  are  constant  or  Gaussian,  and  also  for  some 
other  very  particular  cases,  such  as  for  gamma  random  variables  with  the 
same  scale  parameter.  Otherwise,  it  is  only  an  approximation  that  is  not 
necessarily  accurate  in  all  cases.  For  instance,  Stewart  et  al.  (2007)  ob¬ 
serve  that  this  approach  is  only  reasonably  accurate  when  summing  two 
gamma  random  variables  if  the  shape  parameters  are  not  lower  than  0.1 
and  the  scale  parameters  do  not  differ  by  more  than  a  factor  of  10.  Accura¬ 
cy  of  approximating  the  sum  of  two  lognormal  random  variables  by  anoth¬ 
er  lognormal  random  variable  may  be  adequate  in  one  tail  of  the  summed 
distribution  but  not  the  necessarily  the  other  (Mehta  et  al.  2007).  The  ac¬ 
curacy  of  this  approximation  with  other  distributions  such  as  the  trans¬ 
formed  Rice-Nakagami  is  unverified.  While  the  method  seems  to  make  in¬ 
tuitive  sense,  the  accuracy  of  the  approach  is  difficult  to  confirm  or  predict 
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for  various  distributions  and  cases.  It  is  reasonable  to  assume  that  there 
are  limitations  and  that  it  is  not  good  in  all  cases. 

However,  this  method  is  a  very  simple  alternative  to  other  more  robust  yet 
computationally  intensive  approaches  and  much  easier  to  implement.  It  is 
also  important  to  note  that  parametric  pdfs  representing  various  signals 
are  themselves  only  estimates  to  begin  with  because  both  the  models  and 
the  parameters  in  the  models  used  to  generate  them  are  approximate  due 
to  idealizations  of  the  signal  physics  and  uncertainties  in  the  weather, 
ground  conditions,  source  signature  and  directivity,  noise  background,  etc. 
Especially  if  it  will  greatly  increase  computational  intensity,  it  is  not  neces¬ 
sarily  worthwhile  to  compute  the  sum  of  multiple,  imprecise  random  va¬ 
riables  with  very  high  precision. 

Although  the  exact  sum  of  multiple  random  variables  of  related  type  can¬ 
not  typically  be  represented  by  an  analytical  pdf,  this  method  essentially 
suggests  selecting  one  that  will  serve  as  a  good  approximation  with  mean 
and  variance  conserved.  Generally,  the  sum  is  best  approximated  by  the 
pdf  of  one  of  the  summed  random  variables,  but  sometimes  it  is  better  to 
use  a  pdf  that  is  different  but  still  a  generalized  form  of  the  summed  ran¬ 
dom  variables.  For  example,  the  sum  of  multiple  exponential  random  va¬ 
riables  would  be  better  approximated  by  a  gamma  distribution  rather  than 
another  exponential  distribution,  where  means  and  variances  are  con¬ 
served.  In  fact,  this  would  be  exact  if  the  exponential  random  variables 
were  all  identical  (with  the  same  mean).  Likewise,  the  summation  of  mul¬ 
tiple  random  variables  of  different  but  related  type  is  represented  as  a  pdf 
of  the  most  generalized  form  of  the  related  random  variables.  An  example 
of  this  is  the  sum  of  an  exponential  random  variable  and  a  gamma  random 
variable,  which  would  be  a  gamma  random  variable  with  the  means  and 
variances  conserved.  Here,  the  exponential  and  gamma  distributions  are 
closely  related,  belonging  to  the  exponential  class  of  density  functions, 
where  the  gamma  distribution  with  shape  parameter  equal  to  1  reduces  to 
the  exponential  distribution  with  a  mean  equal  to  the  scale  parameter  (i.e., 
the  exponential  distribution  is  a  special  case  of  the  gamma  distribution). 
So,  in  fact,  the  sum  of  exponential  and  gamma  random  variables  is  essen¬ 
tially  just  a  sum  of  gamma  random  variables. 

Conceivably  (but  not  necessarily),  the  method  could  well  approximate  the 
sum  of  many  different  kinds  of  random  variables,  so  long  as  a  parametric 
distribution  that  generalizes  a  diverse  set  of  random  variables  is  available. 
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For  instance,  the  sum  of  any  combination  of  two  or  more  of  the  lognormal, 
exponential,  Weibull,  and  gamma  distributions  may  be  well  approximated 
by  a  generalized  gamma  distribution  with  means  and  variances  conserved, 
since  the  aforementioned,  well-known  distributions  are  actually  special 
cases  of  the  generalized,  three-parameter  gamma  distribution.  Indeed, 
there  are  other  three-parameter  distributions  that  generalize  various  well- 
known  distributions  such  as  the  generalized  and  skew  normal  distribu¬ 
tions,  both  of  which  can  mimic  other  skewed  distributions  like  the  gamma, 
lognormal,  and  Weibull  distributions  but  also  includes  the  normal  distri¬ 
bution  as  a  special  case. 

Though  there  is  an  added  complication  in  representing  the  sum  of  random 
variables  with  a  three-parameter  distribution,  since  the  three  parameters 
must  be  estimated  numerically  via  maximum  likelihood  or  method  of 
moments,  it  is  possible  and  potentially  worthwhile  if  the  process  proves  to 
be  relatively  less  computationally  intensive  (Stacy  and  Mihram  1965; 
Ashkar  et  al.  1988;  Varanasi  and  Aazhang  1989).  When  using  the  method 
of  moments,  it  is  helpful  to  notice  that  the  third  unstandardized  central 
moment  can  be  computed  exactly  for  the  sum,  since  they  are  also  “con¬ 
served”  like  the  means  and  variances,  which,  in  fact,  are  the  first  raw  mo¬ 
ment  and  the  second  unstandardized  central  moment,  respectively.  This 
results  in  three  moments  that  can  be  computed  exactly  for  the  sum,  which 
can  be  used  to  estimate  the  required  three  parameters  to  obtain  the  gene¬ 
ralized  distribution.  Otherwise,  it  may  be  advantages  to  compute  a  few 
particular  cumulants  of  the  generalized  distribution  of  the  sum  to  estimate 
its  parameters  instead,  since  all  cumulants  are  conserved. 

As  described,  this  method  only  applies  for  summing  multiple  random  va¬ 
riables  of  related  type.  If  no  generalized  distribution  for  two  or  more  di¬ 
verse  random  variables  is  available,  it  is  impossible  to  predict  the  appro¬ 
priate  analytical  pdf  that  best  approximates  their  sum.  For  example,  the 
lognormal  and  the  transformed  Rice-Nakagami  random  variables  are  not 
of  related  type  and  do  not  have  a  parametric  distribution  that  generalizes 
both  of  them,  so  this  method  does  not  apply  in  this  case. 

Finally,  this  method  is  both  commutative  and  associative  in  that  the  order 
in  which  multiple  random  variables  are  summed  together  does  not  affect 
the  final  representation  of  their  sum.  This  follows  first  from  the  fact  that 
the  summation  of  means  and  variances  (and  the  third  unstandardized  cen¬ 
tral  moment,  if  necessary)  is  itself  commutative  and  associative.  Then,  the 
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method  selects  the  analytical  pdf  of  the  same  type  or  of  a  more  generalized 
type  than  that  which  describes  the  two  summed  random  variables,  regard¬ 
less  of  which  one  is  presented  first  in  the  method’s  algorithm. 

Using  the  Laplace  approximation 

The  Laplace  approximation  is  another  approach  that  introduces  the  ap¬ 
proximation  at  a  different  stage  in  the  summation  of  random  variables. 
Rather  than  approximating  the  resulting  sum  of  multiple  random  va¬ 
riables  with  an  appropriate  analytical  pdf  as  in  the  previously  described 
method,  it  is  possible  to  approximate  each  non-Gaussian  random  variable 
addend  by  a  Gaussian  pdf.  Then,  the  sum  of  these  Gaussian  approxima¬ 
tions  may  be  computed  exactly,  since  the  exact  sum  of  multiple  Gaussian 
random  variables  is  simply  another  Gaussian  random  variable  with  means 
and  variances  conserved. 

Since  many  familiar,  unimodal  random  variables  may  be  roughly  approx¬ 
imated  by  a  Gaussian  pdf,  the  Laplace  approximation  is  widely  used  in 
many  situations  requiring  manipulation  of  multiple  random  variables  that 
must  be  or  is  greatly  preferred  to  be  Gaussian.  In  our  case,  obtaining  a 
good  Gaussian  approximation  of  multiple  random  variables  greatly  simpli¬ 
fies  their  summation  by  avoiding  the  use  of  a  more  computationally  inten¬ 
sive,  numerical  approach.  It  is  also  likely  that  the  Laplace  approximation  is 
within  the  precision  errors  inherent  in  the  idealizations  and  uncertainties 
of  the  models  and  parameters  used  in  simulating  the  signal-transmission 
and  sensing  process. 

Since  the  Laplace  method  involves  the  second-degree  Taylor  expansion  of 
the  given,  non-Gaussian  pdfs  logarithm  centered  on  the  mode,  both  the 
mode  and  the  second  derivative  of  the  pdfs  logarithm  at  the  mode  must  be 
computed.  Often  these  can  be  calculated  with  simple,  closed-form  expres¬ 
sions.  If  not,  quick  numerical  computations  can  be  devised  noting  that  the 
mode  is  defined  as  where  the  maximum  of  the  pdf  occurs.  Details  on  the 
theory  and  definition  of  the  Laplace  approximation  can  be  found  in  Bishop 
(2006). 

The  convolution  approach 

Though  it  is  the  most  computationally  intense,  it  is  possible  to  always 
compute  the  pdf  of  the  sum  of  two  independent  random  variables  via  the 
convolution  of  their  pdfs.  Since  all  pdfs  are  compactly  supported  and 
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locally  integrable,  the  convolution  of  two  pdfs, /and  g,  denoted  as/*  g,  is 
the  following  well-defined  and  continuous  integral  transform: 


def  - 

The  preferred  method  for  computing  the  convolution  is  via  the  forward 
and  inverse  discrete  Fourier  transforms,  for  which  fast  Fourier  transform 
algorithms  are  available.  According  to  the  convolution  theorem,  the  Fouri¬ 
er  transform  of  a  convolution  is  the  pointwise  product  of  the  Fourier  trans¬ 
forms  of  the  two  functions  being  convoluted  together: 

F[f*g]  =  k-F[f]-F[g\ 

where  k  is  a  constant  that  depends  on  the  normalization  of  the  Fourier 
transform  and  f\j]  denotes  the  Fourier  transform  of/  Then,  the  inverse 

Fourier  transform  may  be  applied  to  obtain  the  convolution,/*  g: 

f  *  g  =  F  '  W  *  g]]  =  F  1  \k  ■  F\f]  ■  F[g]]  =  kF~'  [F\f]  ■  F\g\ 

Since  the  characteristic  function  of  a  random  variable  is  a  Fourier  trans¬ 
form  of  its  pdf  and  is  available  for  every  type  of  random  variable,  it  is  then 
possible  to  compute  the  sum  of  multiple  random  variables  using  the 
inverse  Fourier  transform  and  characteristic  functions  by  the  general 
formula  above. 

By  using  the  convolution  theorem  and  the  fast  Fourier  transform  (FFT) 
algorithms,  an  implementation  for  computing  the  sum  of  multiple  random 
variables  via  the  numerical  convolution  of  two  pdfs  may  be  made  practical. 
In  practice,  the  FFT  algorithm  must  be  implemented  carefully  to  avoid 
errors  from  discretization  and  truncation  of  integrals  as  much  as  possible. 
This  can  be  quite  difficult,  but,  if  done  correctly,  the  convolution  approach 
would  make  it  possible  to  sum  multiple  random  variables  of  even  vastly 
different  type.  It  is  also  commutative  and  associative,  as  required. 

Central  limit  theorem 

As  more  independent  random  variables  are  summed,  the  joint  random 
variable  monotonically  converges  into  a  Gaussian  distribution  given 
certain  conditions.  By  the  classical  central  limit  theorem,  the  sum  must  be 
of  many  independent  and  identically  distributed  random  variables,  each 
having  positive  variance  with  both  mean  and  variance  being  finite. 
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However,  the  theorem  actually  also  holds  for  non-identically  distributed 
random  variables  if  Lyapunov’s  condition  is  satisfied,  which  further 
implies  that  Lindeberg’s  condition  is  also  met  (Ash  and  Doleans-Dade 
2000).  Though  the  rate  of  convergence  depends  on  the  type(s)  of  the 
identical  or  non-identically  distributed  random  variable  addends,  it  is 
conceivable  that  the  sum  of  even  three  to  five  non-Gaussian  random 
variables  may  already  be  well  approximated  by  a  Gaussian  distribution. 
Hence,  in  cases  where  enough  random  variables  are  added  together,  where 
at  least  one  of  them  is  non-Gaussian,  the  sum  may  be  simply  represented 
by  a  Gaussian  distribution  with  means  and  variances  conserved.  Any 
implementation  for  summing  multiple  random  variables  should  take  this 
limit  theorem  into  account  to  improve  both  accuracy  as  well  as  com¬ 
putational  speed.  When  the  theorem  applies,  the  sum  should  be  directly 
represented  as  Gaussian  with  means  and  variances  conserved  rather  than 
computing  many  Laplace  approximations,  numerically  evaluating  many 
convolutions,  or  using  a  different  pdf  for  the  final  representation. 
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5  Numerical  Issues 

Although  most  statistical  operations  required  for  the  pdfs  described  in  this 
report  are  relatively  straightforward  to  program,  it  is  helpful  to  illustrate 
some  of  the  operations  that  make  use  of  certain  numerical  algorithms.  It  is 
also  worthwhile  to  describe  how  more  complicated  operations  must  be 
programmed  carefully  so  that  they  compute  properly  using  floating-point 
arithmetic  on  a  computer.  For  the  latter  concern,  only  a  selection  of  exam¬ 
ples  is  presented  to  illustrate  some  common  problems  that  may  arise.  The 
Appendix  includes  further  details. 

Computing  the  quantile  using  the  bisection  method 

With  the  exception  of  the  exponential  model,  for  which  a  closed-form  for¬ 
mula  for  evaluating  the  quantile  is  available,  all  other  statistical  models  are 
coded  to  numerically  solve  for  the  quantile,  x,  by  approximating  the  root, 
y0,  of  the  following  function,  Q(y) ,  using  the  bisection  method  for  a  given 
p  =  cdf(x)  with  the  tolerance  set  to  the  lowest  representable  number  in  bi¬ 
nary  double  (64-bit)  floating-point  precision  (approximately -1.8  x  io3°8): 


Q(y)  =  cdf  {y)-P 

where  the  root,  y0,  is  equal  to  the  quantile,  x: 

Q(}’o )  =  cdf(v0  )-P  =  cdf  (y0 )-  cdf  (x)  =  0  <=>  cdf  (y0 )  =  cdf  (x)  <=>  y0  =  x 

The  lower-  and  upper-bounds  used  for  the  bisection  method  is  set  to  span 
the  entire  domain  supported  by  a  particular  statistical  distribution  that  is 
representable  in  binary  double  (64-bit)  floating-point  precision.  Specifical¬ 
ly,  the  lower-  and  upper-bounds  for  the  lognormal,  gamma,  and  x  ->  x2 
transformed  Rice-Nakagami  models  are  set  to  zero  and  the  highest  repre¬ 
sentable  number  (approximately  1.8  x  io3°8).  For  random  variables  with 
real-number  support  like  the  Gaussian,  the  lower-  and  upper-bound  is  set 
to  the  lowest  and  highest  representable  number  (approximately  -1.8  x 
103°8  and  +1.8  x  io3°8),  respectively.  The  bisection  method  used  to  solve 
for  the  quantile  requires  that  the  cdf  be  always  computable  within  the 
bounds.  Since  the  bounds  are  set  to  span  very  large  domains,  the  cdf  must 
be  programmed  carefully  so  that  it  remains  computable  for  all  input  values 
within  the  domain.  Concerns  related  to  how  the  cdf  should  be  pro¬ 
grammed  to  satisfy  this  required  condition  are  discussed  in  the  last  two 
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sections  of  this  chapter  (“Numerical  stability  and  well-posed  problems” 
and  “Illustrative  examples  of  numerically  stable  statistical  computations”). 

Although  the  bisection  method  only  converges  linearly  (i.e.,  the  absolute 
error  of  the  true  and  estimated  roots  is  halved  at  each  step),  it  is  preferred 
here  because  of  it  robustness.  Unfortunately,  other  methods  with  higher 
convergence  rates  such  as  the  Newton  and  Secant  methods  fail  due  to  the 
nature  of  the  cdf  curve,  which  is  asymptotically  a  horizontal  line  (with  the 
first  derivative  equal  to  zero)  at  both  ends.  An  optimized  algorithm  with  a 
mixture  of  different  methods  is  conceivable  but  does  not  necessarily  im¬ 
prove  the  convergence  rate  appreciably. 

Computing  the  cdf  using  Gauss-Legendre  quadrature  on  the  pdf 

Usually  there  is  an  expression  for  the  cdf  that  could  be  directly  imple¬ 
mented,  but,  if  it  is  complicated  (e.g.,  it  contains  a  special  function  that  is 
difficult  to  program),  numerical  integration  of  the  pdf  may  be  a  simpler 
alternative,  although  it  is  generally  more  computationally  intensive: 


x 

cdf(x)  =  J  pdf(/)# 

0 

The  Gauss-Legendre  quadrature  rule  is  a  very  efficient  method  for  per¬ 
forming  the  integration.  Basically,  the  method  approximates  the  definite 
integral  of  a  function  via  a  weighted  sum  of  function  values  at  specified 
points  within  the  domain  of  integration.  More  details  about  stably  imple¬ 
menting  this  method  for  arbitrary  integration  domains  are  provided  in  the 
Appendix. 

In  addition  to  the  lower-  and  upper-bounds  of  the  integration,  the  method 
requires  specification  of  the  number  of  weighted  function  values,  N,  to  be 
summed.  The  appropriate  number  of  points  for  the  Gauss-Legendre  rule, 
N,  so  that  the  integration  of  the  pdf  is  computed  to  a  sufficient  level  of  pre¬ 
cision  depends  on  the  nature  of  the  pdf  curve  over  the  specified  range  of 
integration.  Since  the  theory  behind  the  AT-point  Gauss-Legendre  rule  en¬ 
sures  that  it  computes  exact  integrals  for  the  class  of  polynomials  of  de¬ 
gree  2N-1  or  less,  one  must  essentially  predict  the  degree  of  the  polynomial 
that  adequately  approximates  the  pdf  over  the  specified  range  of  integra¬ 
tion  to  find  the  appropriate  N.  As  N  is  increased,  the  numerical  integration 
becomes  more  computationally  intensive  (and  time-consuming),  so  it  is 
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important  to  choose  an  N  that  is  not  greater  than  what  suffices  to  obtain 
results  of  sufficient  precision. 

To  minimize  the  number  of  points  necessary  in  the  Gauss-Legendre  rule 

x 

and  still  compute  cdf  (x) = |  pdf(r  with  sufficient  precision,  it  is  helpful  to 

0 

minimize  the  range  of  integration  as  much  as  possible  by  programming  cdf 
computation  differently  depending  on  x: 


cdf(x)  =  < 


0,  for  pdf  (x)  =  0  and  x  <  mean 

X 

| pdf  (t)dt,  for  pdf(x)  >  0  and  x  <  mean 

/ 


u 

1  - 1  pdf(/)rff,  for  pdf  (x)  >  0  and  x  >  mean 

x 

1,  for  pdf(x)  =  0and  x  >  mean 


where  /  and  u  are  approximations  of  the  lowest  and  highest  values  for  x 
with  nonzero  pdf  values,  respectively.  More  precisely,  l  and  u  are  com¬ 
puted  as  the  first  integer  standard  deviate  below  and  above  the  mean  with 
a  pdf  value  of  zero,  respectively.  If  /  and/or  u  is  beyond  the  given  statistical 
distribution’s  representable  support  domain,  it  is  set  as  the  lowest  and 
highest  representable  value  of  the  support  domain.  Then,  empirical  inves¬ 
tigations  suggest  that  the  50-point  Gauss-Legendre  rule  would  be  suffi¬ 
ciently  precise  for  computing  any  desired  cdf  value  (Figure  3).  Although  all 
parameters  cannot  be  tested,  it  is  reasonable  to  assume  that  the  general 
behavior  of  all  pdfs  of  the  same  type  is  similar  enough  that  the  50-point 
rule  would  always  suffice. 


This  approach  is  currently  used  for  the  X  -» X1  transformed  Rice- 
Nakagami  model  to  avoid  implementing  the  first-order  Marcum  Q- 
Function  in  its  cdf  expression.  It  is  also  included  in  the  gamma  model  as 
an  alternative  in  the  potential  event  that  the  implementation  for  the  direct 
approximation  fails  for  certain  inputs.  As  it  turns  out,  some  care  must  be 
taken  in  programming  floating-point  implementations  of  certain  pdfs  and 
cdfs  (including  the  gamma  cdf),  as  will  be  discussed  in  the  following  sec¬ 
tions. 
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Figure  3.  Comparison  of  cdf  generated  using  different  Appoint  Gauss-Legendre  rules  to 
integrate  the  A"— ►  A2  transformed  Rice-Nakagami  pdf  with  mean  =  100  and  var  =  70. 

Vertical  axis:  pdf(x);  horizontal  axis:  x.  From  left  to  right  on  the  top  row:  N  =  2,  5,  and  10.  From 
left  to  right  on  the  bottom  row:  N  =  20,  50,  and  100.  As  N  is  increased,  accuracy  of  the  pdf 
integration  is  improved  and  a  smooth  cdf  curve  is  generated. 

Numerical  stability  and  well-posed  problems 

When  programming  the  various  statistical  operations,  it  is  important  to 
code  calculations  so  that  they  function  properly  when  performed  using 
floating-point  arithmetic  on  a  computer.  For  programming  very  simple 
calculations,  it  is  not  necessarily  important  to  understand  the  mechanics 
of  floating-point  arithmetic.  It  is,  however,  essential  to  understand  when 
formulating  more  complicated  functions  and  numerical  algorithms.  When 
a  calculation  is  programmed  to  function  properly  within  the  limitations  of 
floating  point,  it  is  called  numerically  stable.  Goldberg  (1991)  has  a  much 
more  comprehensive  treatment  of  the  subject  than  that  which  is  provided 
here. 

The  basic  problem  is  caused  by  computers  with  limited  memory  that  can 
only  represent  (and,  therefore,  work  with)  a  finite,  discrete  set  of  numbers. 
Floating  point  describes  a  system  for  allocating  computer  memory  effi¬ 
ciently  so  that  a  very  large  set  of  numbers  can  be  made  “representable.”  A 
more  thorough  background  on  floating  point  is  provided  in  the  Appendix. 

A  variety  of  issues  due  to  floating  point  are  well  known.  First,  there  is 
round-off  error,  which  is  the  difference  between  the  calculated  approxima¬ 
tion  of  a  number  and  its  exact  mathematical  value.  This  type  of  error  is  in- 
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troduced  when  the  inputs  and  outputs  of  an  arithmetic  operation  must  be 
rounded  to  the  nearest  representable  number.  Typically,  round-off  errors 
are  insignificant  in  simple  computations,  except  in  relatively  rare  instances 
(e.g.,  when  two  large,  nearly  equal  numbers  are  subtracted,  causing  the 
number  of  accurate  significant  digits  in  the  difference  to  be  greatly  re¬ 
duced,  sometimes  called  a  catastrophic  cancellation ).  But,  small  round-off 
errors  may  accumulate  unacceptably  when  performing  operations  that  re¬ 
quire  a  sequence  of  many  calculations  (e.g.,  large  summations,  matrix  in¬ 
version,  eigenvector  computation,  and  solving  differential  equations). 

Another  concern  is  more  critical  in  our  case,  particularly  when  program¬ 
ming  the  pdf,  cdf,  and  quantile.  Given  that  the  final  result  of  a  computa¬ 
tion  is  indeed  representable,  it  is  important  that  the  results  of  the  various 
arithmetic  operations  within  that  computation  remain  within  the  given 
floating-point  representation’s  domain  of  representable  numbers.  Ideally, 
a  computation  should  be  programmed  so  that  this  condition  is  met  for  all 
representable  inputs.  When  this  condition  is  not  met,  the  final  result  is 
typically  either  very  inaccurate  or  not  computable  at  all,  since  an  arithmet¬ 
ic  operation  has  resulted  in  an  overflow  or  underflow  (which  is  denoted  by 
saying  that  the  computation  has  suffered  from  premature  overflow  or  un¬ 
derflow,  respectively).  Then,  the  computation  could  unpredictably  return 
an  undefined  or  infinite  number  (e.g.,  +oo,  -oo,  or  NaN),  although  the  true, 
mathematical  solution  to  the  computation  is  within  the  representable 
range.  For  example,  computation  of  the  average  of  two  nonzero  numbers, 
x  and  y,  that  are  within  the  representable  domain  using  the  formula, 

(x  +  y)  /  2,  may  return  +oo  or  -oo  if  (x  +  y)  overflows,  although  the  average 
would  always  lie  within  the  representable  range. 

To  address  this  latter  issue,  it  is  helpful  to  notice  that  sometimes  a  single 
calculation  can  be  achieved  in  several  ways.  In  the  previous  example  for 
computing  the  average  of  the  two  numbers,  x  and  y,  it  is  better  to  program 
the  distributed  formula,  x  /  2  +  y  /  2,  which  does  not  have  the  overflow 
problem.  But,  this  formula  is  now  susceptible  to  underflow  by  the  terms, 
x  /  2  or  y  /  2,  for  small  enough  x  or  y.  Since  premature  underflow  would 
only  present,  at  most,  an  absolute  error  of  twice  the  smallest  representable 
number,  it  is  more  important  to  program  against  premature  overflow  than 
premature  underflow  in  this  case— for  other  situations  when  underflow 
may  present  a  division-by-zero  operation,  premature  underflow  is  much 
more  important  to  prevent.  The  final,  numerically  stable  algorithm  for 
averaging  the  two  numbers,  x  and  y,  would  compute  the  distributed 


ERDC/CRREL  TR-10-12 


22 


formula,  x  /  2  +  y  /  2,  only  if  the  original  formula,  (x  +  y)  /  2,  results  in  an 
overflow. 

If  a  computation  is  rather  complicated,  a  careful,  numerically  stable  float¬ 
ing-point  computer  implementation  of  it  can  be  quite  difficult  to  program. 
Although  various  algebraic  manipulations  are  commonly  used,  each  com¬ 
putation  is  unique  and  requires  different  techniques.  While  it  is  easier  to 
program  a  computation  that  must  be  stable  only  within  a  limited  domain 
of  inputs,  the  inputs  that  the  computation  must  be  able  to  handle  may  be 
difficult  to  know  or  predict. 

Illustrative  examples  of  numerically  stable  statistical  computations 

Since  they  share  some  general  similarities,  many  statistical  formulae  suffer 
from  common  issues.  For  instance,  many  pdfs  are  of  the  following  form: 

pdf  (x)  =  ab  exp(-  erf)  =  — ^ — r  • 
exp(crf) 

A  simple,  direct  implementation  of  the  pdf  expression  may  not  normally 
present  issues  with  ordinary  values  for  a,  b,  c,  and  d.  However,  this  formu¬ 
lation  is  certainly  not  robust  for  all  representable  a,  b,  c,  and  d.  Namely, 
the  numerator  and  denominator  may  both  overflow,  although  their  ratio  is 
a  perfectly  ordinary  value  (that  always  lies  within  the  representable  range). 

This  potential  problem  is  easily  remedied  by  a  single,  straightforward  al¬ 
gebraic  rearrangement.  Specifically,  the  division  of  two  potentially  large 
numbers,  A  and  B,  is  better  encoded  via  the  difference  of  their  logarithms: 

^=exp[]n(/f)-]n(/?)] 

D 

so  that  the  ratio  becomes  computable  for  a  much  more  expansive  domain 
of  A  and  B,  since  it  is  only  required  that  their  logarithms  not  overflow.  Us¬ 
ing  this  rearrangement,  the  general  pdf  expression  above  becomes: 


pdf(x)  =  - 


ab 


exp(crf) 


=  exp 


In 


ab 


exp  (erf) 


=  exp[ln(erf>)  -  ln(exp(crf))]  =  exp[ln(aZ>)  -  erf]  =  exp[ln(a)  +  ln(&)  -  erf] 


It  is  important  to  note  here  that  not  all  overflows  are  necessarily  an  issue 
but  that  they  become  problematic  only  when  inaccuracies  may  result.  For 
example,  the  potential  overflow  of  the  term,  cd,  in  the  above  rearrange¬ 
ment  is  not  an  issue  because  it  would  always  imply  that  the  pdf  is  zero 
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and,  in  fact,  actually  results  in  a  zero  to  be  returned  for  the  pdf.  (It  is 
helpful  to  notice  here  that  the  rearranged  pdf  would  always  equal  zero 
when  the  term,  cd,  becomes  large  enough  even  before  it  overflows.)  Over¬ 
flow  of  A  and/or  B  in  the  previous  expression  is  troublesome  because  it 
does  not  imply  anything  about  the  true  value  of  A  /  B  and  inaccurately 
returns  o,  +oo,  or  NaN. 


Since  computing  exponentials  and  logarithms  are  more  computationally 
expensive,  it  is  better  to  calculate  the  stable  rearrangement  as  an  alterna¬ 
tive  only  when  the  floating-point  computation  of  A  /  B  fails,  if  the  simpler 
formulation  is  only  rarely  unstable.  In  certain  cases,  however,  it  may  be 
that  the  straightforward  formula,  A  /  B,  is  actually  more  often  unstable 
than  it  is  stable.  This  is  true,  for  example,  with  the  gamma  cdf,  which  will 
be  described  shortly.  In  these  cases,  it  would  be  more  efficient  to  just  com¬ 
pute  the  stable  rearrangement  directly.  The  Appendix  includes  specific  de¬ 
tails  on  how  to  apply  this  rearrangement  for  each  type  of  distribution. 


As  just  stated,  this  technique  is  useful  in  coding  the  gamma  cdf,  defined 
for  a  >  o ,  a  >  o,  and  (3  >  o: 


f 

cdf(x)  =  P 

V 


r(«) 


where  P  is  the  regularized  incomplete  gamma  function: 


n(  \  AS’X) 

p\s,x)=fn-r- 

rb) 


1  x 

1  f  s-\ 

rb){ 


'  dr 


and  T  and  y  are  the  gamma  function  and  the  lower  incomplete  gamma 
function,  respectively. 


In  this  formula,  the  individual  terms, 


a 1  and  r(«),  overflow  at  modest 

V  P) 


values  of  x,  a,  and  f>,  although  the  regularized  incomplete  gamma  function 
itself  has  the  limiting  values: 


^b,o)  =  o  and  p(?,°o)=i. 


Thus,  it  is  also  helpful  in  this  case  to  implement  the  exponentiation  of  the 
difference  of  the  logarithms  of  the  two  terms: 
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(  \  A a' 

cdf(x)  =  P  a,—  = - T 

l  P)  r(£ 


P 

a) 


=  exp 


r 


In 


a’  P 

r(«) 


:  exp 


In 


(  f 


a, 

v  Pjj 


-In(r(«)) 


Though  they  are  used  less  frequently,  various  other  manipulations  are  use¬ 
ful  in  more  specific  circumstances  (and  often  in  combination  with  the  first 
technique).  For  instance,  after  using  this  first  technique,  the  Gaussian  cdf 
still  contains  a  term  that  resembles  the  following: 


ln(a  +  b )  • 


For  all  representable  a  and  b,  the  true  value  of  this  term  is  always  within 
the  representable  range,  but  +oo  would  be  returned  whenever  (a  +  5)  over¬ 
flows.  This  issue  can  be  avoided  by  programming  the  following: 


ln(a  +fe)  =  ln 


where  c  is  set  as  the  largest  representable  number  in  the  given  floating¬ 
point  system. 


A  similar  distributive  approach  is  used  in  multiple  other  situations,  includ¬ 
ing  the  formula  for  the  gamma  pdf,  which  presents  the  following  general 
expression: 

ab  —  c . 

When  implemented  directly,  the  subtraction  by  c  is  ignored  whenever  the 
term,  ab,  overflows  and  results  in  the  entire  expression  being  automatical¬ 
ly  equated  to  +oo  or  -oo.  This  is  unfortunate  if  c  is  of  the  same  sign  as  ab, 
since  the  entire  expression  may  very  well  be  within  the  representable 
range  (if  c  is  of  large  enough  magnitude).  This  problem  can  be  easily 
avoided  by  programming  the  following  instead: 


ab-c  =  d\ 


ab  -  c) 
d  ) 


where  d  is  selected  to  be  large  enough  so  that  the  term,  — ,  does  not  over¬ 
flow. 

In  many  instances,  it  is  only  important  to  simply  pay  attention  to  the  order 
of  operations.  This,  in  fact,  applies  for  programming  the  term,  ,  in  the 
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technique  explained  above.  It  also  applies  in  many  other  situations,  in¬ 
cluding  when  programming  a  term  in  the  lognormal  cdf,  which  is  of  the 
following  form: 

a 

be 

For  even  a  simple  expression  like  this,  there  are  many  ways  to  program  it: 


a  _  a 
be  (be) 


With  some  quick  investigation,  it  is  relatively  easy  to  confirm  that  (usually) 
only  one  of  these  implementations  (if  any)  is  stable  given  certain  limita¬ 
tions  on  the  domain  of  a,  b,  and  c.  But,  in  a  couple  of  special  cases  of  this 
general  expression,  it  is  always  best  to  implement  them  in  the  following 
way: 

a  _\b 
b2  b 


and 


Simply  changing  the  order  of  operations  is  often  sufficient  throughout 
parts  of  various  statistical  formulae;  however,  attention  to  this  detail  is 
quite  easily  overlooked.  Application  of  this  technique  in  specific  cases  is 
described  within  the  Appendix. 


When  no  particular  ordering  of  operations  is  universally  stable  for  the  giv¬ 
en  domains  of  a,  b,  and  c,  it  is  always  possible  to  resort  to  the  following 
general  rearrangement  that  is  stable  for  all  representable  a,  b,  and  c: 


a 

—  =  exp 
be 


=  exp[ln(a)  -  ln(i>)  -  ln(c)]  • 


Finally,  it  is  helpful  to  note  certain  consequences  of  finite-precision  arith¬ 
metic  when  analyzing  and  implementing  stable  mathematical  expressions 
in  floating  point.  Essentially,  any  floating-point  system  represents  a  finite, 
discrete  set  of  numbers  along  the  real  number  line  that  retains  the  same 
relative  difference  only  between  each  representable  number  (Figure  4). 
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This  results  in  additions  or  subtractions  by  sufficiently  small  numbers  to 
be  ignored  because  the  floating-point  system  is  not  precise  enough  to 
represent  the  difference.  This  observation  is  useful  when  programming, 
for  instance,  the  following  general  expression,  which  arises  in  lognormal 
parameter  conversions: 

ln(a  +  be) 

which  becomes  the  following  for  sufficiently  large  differences  in  magni¬ 
tude  between  a  and  be: 


ln(a  +  be) — >  In  (be)  =  ln(/>)  +  ln(c) 


This  is  a  simpler  alternative  to  the  earlier,  distributive  technique: 


ln(a  +  b)  =  In 


if  it  can  be  known  that  there  will  always  be  a  sufficiently  large  magnitude 
difference  between  the  added  or  subtracted  terms. 


This  observation  is  often  useful  in  analyzing  the  stability  of  an  implemen¬ 
tation  and  concluding  that  certain  potential  underflow  issues  are  accepta¬ 
ble  to  ignore.  Specifically,  the  magnitude  difference  of  a  given  representa¬ 
ble  number  and  an  underflowing  number  may  be  large  enough  that  the 
underflow  to  zero  would  not  affect  their  summation,  since  the  contribution 
by  the  underflowing  number  would  be  ignored  anyways  even  if  it  was  not 
necessarily  underflowing.  However,  there  are  situations  when  it  is  prudent 

r —  denormal — , 


under-  under¬ 
flow  flow 

overflow  usable  range  usable  range  overflow 
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to  check  and  fix  a  complicated  implementation  for  potential  underflow 
issues,  particularly  if  the  underflow  would  result  in  a  division  or  multipli¬ 
cation  by  zero.  For  example,  it  would  be  better  to  program  the  following 
two  expressions  as  shown: 


( ab\cd )  =  exp[ln((a6)(a/ ))]  =  exp[ln(a&)  +  ln(crf )]  =  exp[ln(a)  +  ln(/>)  +  ln(c)  +  ln(t/ )] 


and 


M 

{cd) 


=  exp 


In 


{cd  ) 


=  exp[ln(aZ?)  -  ln(cc/)]  =  exp[ln(a)  +  ln(/>)  -  ln(c)  -  ln(rf)] 


so  that  underflow  by  the  terms,  ab  and/or  cd,  would  not  result  in  potential 
inaccuracies  due  to  the  expressions  being  automatically  computed  as  o, 
+oo,  -oo,  or  NaN.  The  Appendix  includes  more  detailed  explanations  of  po¬ 
tential  underflow  issues  for  particular  statistical  computations. 


There  are  obviously  many  more  algebraic  manipulations  that  could  be  use¬ 
ful  in  different  situations;  the  listed  techniques  here  are  by  no  means 
exhaustive.  They  are,  however,  what  is  most  helpful  for  carefully  pro¬ 
gramming  various  statistical  operations  in  floating  point.  As  the  mechan¬ 
ics  of  floating  point  are  better  understood,  it  becomes  easier  to  discover 
other  appropriate  techniques  and  improve  coding  habits.  Gradually,  it 
becomes  obvious  that  certain  expressions  should  be  coded  always  in  a 
certain  way  such  as: 


ln(aZ?)  =  ln(a)  +  ln(Z?) 

and 

sfab  =  4a4b 

As  it  turns  out,  designing  stable  floating-point  implementations  is  highly 
nontrivial.  But,  it  is  essential  if  the  domain  of  inputs  is  expansive  or 
unknown  and/or  if  a  particular  numerical  algorithm  depends  on  it  to  func¬ 
tion  properly.  In  the  case  here,  for  example,  stable  implementations  of  the 
cdf  is  vital  if  the  algorithm  that  solves  for  the  quantile  uses  a  root-finding 
method  with  bounds  set  to  span  a  very  expansive  domain  (e.g.,  the  entire 
representable  region).  Only  a  broad,  conceptual  discussion  is  given  here 
for  brevity.  The  Appendix  gives  specific  details  for  implementing  a  particu¬ 
lar  computation. 
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6  Conclusions 

Given  the  randomness  of  most  signal-generation  and  propagation 
processes,  predictions  of  sensor  performance  are  inherently  probabilistic 
in  nature.  While  it  is  important  to  minimize  (e.g.,  using  accurate  parame¬ 
ters  and  models),  various  types  and  levels  of  uncertainty  would  still  exist 
in  virtually  all  situations  when  modeling  a  complicated  process  such  as 
signal  transmission  and  sensing.  Especially  when  small  variations  in  the 
environmental  state  dramatically  alter  outcomes,  calculating  and  commu¬ 
nicating  predictions  in  a  statistical  manner  is  critical. 

This  report  focused  on  representing  target  and  noise  signals  as  random 
variables  to  account  for  variations  in  sensor  data  (e.g.,  sound  level,  seismic 
energy,  concentration  of  a  chemical  agent,  etc.)  that  are  essentially  sto¬ 
chastic  due  to  many  unknown  and  potentially  dynamic  environmental  var¬ 
iations  affecting  signal  generation  and  propagation.  By  modeling  signals 
statistically,  the  report  discussed  howto  make  probabilistic  predictions  on 
sensor  performance.  Both  the  software  implementation  of  the  statistical 
framework  and  some  important  numerical  issues  in  programming  statis¬ 
tical  computations  were  also  explained. 
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Appendix:  Floating-Point  Implementations  of 
Statistical  Formulae 

This  appendix  includes  specific  details  on  floating-point  implementations 
of  the  various  statistical  formulae  described  in  the  report.  Mathematical 
expressions  are  programmed  so  that  they  always  return  a  representable 
number  whenever  the  true  solution  is  representable  for  all  representable 
parameter  inputs.  When  this  condition  is  met,  it  is  said  that  the  implemen¬ 
tation  is  numerically  stable.  The  representable  numbers  are  often  defined 
by  the  binary  double  (64-bit)  floating-point  precision  system,  which  has 
approximately  16  decimal  digits  of  precision  and  spans  numbers  from 
approximating  -1.8  x  io3°8  to  +1.8  x  io3°8.  A  short  background  on  floating 
point  is  given  as  an  introduction. 

Floating  point 

To  be  able  to  represent  the  widest  range  of  numbers  possible  using  a  fixed 
number  of  bytes  in  computer  memory,  modern  machines  approximate 
numbers  to  a  fixed  number  of  significant  digits  (called  the  significand), 
which  is  scaled  by  multiplication  with  a  base  (normally  2, 10,  or  16)  raised 
to  a  specified  exponent: 


significant  digits  x  baseexpolle”' 

This  numeral  system,  called  floating  point,  allows  the  radix  point  that 
separates  the  integer  part  of  a  number  (to  the  left  of  the  point)  from  the 
fractional  part  (to  the  right  of  the  point)  to  be  placed  anywhere  relative  to 
the  significant  digits  of  the  number  (i.e.,  to  “float”),  making  it  conceptually 
similar  to  scientific  notation.  While  several  different  floating-point  imple¬ 
mentations  have  been  used  in  the  past,  the  Institute  of  Electrical  and  Elec¬ 
tronics  Engineers  (IEEE)  has  standardized  the  practice  in  IEEE  754,  which 
is  now  followed  by  almost  all  modern  computers. 

The  IEEE  754  Standard  for  Floating-Point  Arithmetic  defines  several  basic 
floating-point  formats  using  different  radixes  (i.e.,  bases)  and  amounts  of 
computer  bits,  of  which  two  are  most  widely  used  in  computer  hardware 
and  programming  languages: 


ERDC/CRREL  TR-10-12 


32 


1.  binary  (i.e.,  base  2)  single  precision,  which  occupies  a  total  32  bits  (4  bytes) 
and  has  a  significand  precise  to  24  bits  (about  7  decimal  digits) 

2.  binary  double  precision,  which  occupies  a  total  of  64  bits  (8  bytes)  and  has 
a  significand  precise  to  53  bits  (about  16  decimal  digits). 

Other  basic  floating-point  formats  include  binary  quadruple  precision 
(occupying  128  bits)  and  decimal  (i.e.,  base  10)  formats  encoded  using  64 
or  128  bits.  All  basic  floating-point  formats  also  define  representations  of 
special  values,  including  positive  and  negative  infinities  (+00  and  -00), 
negative  zero  (-0)  distinct  from  ordinary  (“positive”)  zero,  and  a  “not  a 
number”  value  (NaN),  where  NaN  is  generated  by  three  kinds  of  opera¬ 
tions: 

1.  operations  that  produce  indeterminate  forms : 

a.  the  divisions  0/0,  00/00, 00/00,  -00/00,  and  -00/-00 

b.  the  multiplications  ox  00  and  oxoo 

c.  the  power  100 

d.  the  additions  oo+(oc),  (-oo)+oo,  and  equivalent  subtractions  (i.e.,  00-00 
and  -oo-(-oo)) 

2.  real  operations  with  complex  results: 

a.  the  square  root  of  a  negative  number 

b.  the  logarithm  of  a  negative  number 

c.  the  inverse  sine  or  cosine  of  a  number  which  is  less  than  -1  or  greater 
than  +1 

3.  any  operation  with  a  NaN  as  at  least  one  operand 

The  standard  further  defines  howto  simulate  arithmetic  operations  (e.g., 
add,  subtract,  multiply,  divide,  square  root,  exponentiation,  sine,  cosine, 
etc.)  on  floating-point  numbers,  recognizing  certain  exceptional  situations, 
which  include: 

1.  invalid  operations  (e.g.,  square  root  of  a  negative  number) 

2.  division  by  zero 

3.  overflow  (a  number  that  is  too  large  to  represent  correctly) 

4.  underflow  (a  number  that  is  very  small  and  is  inexact) 

5.  inexact  (an  approximated  number) 

Finally,  the  standard  also  defines  various  algorithms  for  rounding  floating¬ 
point  numbers  during  arithmetic  and  conversions,  where  the  default  me¬ 
thod  rounds  to  the  nearest  representable  value  or  to  the  nearest  represent¬ 
able  even  number  if  the  number  lies  midway  between  two  representable 
values. 
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Stable  statistical  implementations 

Details  given  here  for  the  Gaussian,  lognormal  exponential,  and  gamma 
distributions  are  thorough  enough  that  a  reader  may  understand  exactly 
how  they  were  implemented.  Although  the  details  for  the  x  -» x2  trans¬ 
formed  Rice-Nakagami  model  are  not  as  complete,  the  ideas  and  tech¬ 
niques  used  in  its  implementation  are  very  similar  to  those  of  the  other 
models.  The  derivation  for  the  random  variable  transformation  to  obtain 
the  X  -» x1  transformed  Rice-Nakagami  distribution  is  also  provided. 

When  variance  is  set  to  zero,  the  density  function  does  not  exist  and  the 
probability  distribution  becomes  degenerate,  describing  a  discrete  random 
variable  whose  support  consists  of  only  one  value.  While  such  a  distribu¬ 
tion  is  not  random  in  the  practical  sense,  it  does  still  satisfy  the  definition 
of  a  random  variable  and  may  be  described  as  such  in  order  to  provide  a 
way  to  deal  with  constant  values  in  a  probabilistic  framework.  Namely,  the 
probability  mass  function  (pmf)  for  this  case  is  given  by: 

,  .  fl,  for  x  =  mean 
[  0,  otherwise 

Then,  the  cdf  of  the  degenerate  distribution  is  the  Heaviside  step  function 
about  the  mean: 


,  N  fl,  for  x>  mean 

F(x)  =  ' 

[  0,  otherwise 

Here,  “mean”  is  denoted  as  the  single,  constant  value  represented  by  the 
distribution. 

The  described  implementations  in  the  following  sections  do  not  apply 
when  variance  is  set  to  zero.  An  implementation  of  the  probabilistic  de¬ 
scription  explained  above  is  invoked  instead  in  this  case. 

Gaussian  distribution 

Computation  of  the  pdf 

The  following  straightforward  formulation  of  the  Gaussian  pdf ,  defined 
for  x,ju  e  91  and  a2  >0,  is  vulnerable  to  a  variety  of  premature  overflow  and 
underflow  issues,  including  premature  overflow  by  the  exponent  term 
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2cr 


,  which  would, 


in  turn,  cause  premature  underflow  of  the  exponen- 


/ 

tiated  expression,  exp 

V 


(x-v)2) 


where  /u  and  a  are  the  mean  and  standard  deviation  of  the  Gaussian  distri¬ 
bution,  respectively. 

A  simple  reformulation,  however,  eliminates  the  risk  of  premature  over¬ 
flow  and  underflow: 


pdf(x)  =  - 


r  hir 


exp 


pdf(x)  =  exp(ln(pdf  (x)))  =  exp 


v 


rVAr 


exp 


(x-juf 

2cr2 


=  exp 


x  —  JU 


-  ln(var  )+ln(2;r) 


var 


where  the  term  (x  -  /u\  - — —  overflows  only  when  the  true  value  of  the  pdf 


is  zero,  even  when  premature  overflow  occurs  by  the  terms  (x  -  /u)  or 


x  -  jd 
var 


Underflow  of  the  term  ( x-/j 


x-/U 
var  J 


is  not  an  issue  since: 


exp 


(a±min) 


=  1  for  a  6  {F  :  a  *  a  +  min} 


where  F  is  the  set  of  all  numbers  representable  in  binary  double  (64-bit) 
floating-point  precision  and  a  is  set  as  the  expression  ln(var)+ln(27t). 


Computation  of  the  cdf 


Where  the  error  function  is  numerically  stable  for  any  real  number  input 
and  has  the  range  [-1,1],  the  following  formula  for  the  Gaussian  cdf  is  ac¬ 


curate  only  when  computation  of  the  term  -  A  is  numerically  stable: 

(TV  2 


cdf(x)  = 


2 


1  +  erf 


x  —  fi 


A 


7 
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A  numerically  stable  calculation  of  the  term  £  is  done  first  by  changing 

(TV  2 

the  order  of  operations: 


X  ~  jU 

<tV2 


(x~m) 
4i  . 


When  the  term  x-  /u  overflows,  the  term 


x- n 
<tV 2 


is  computed  as  follows: 


X- fl 

aj 2 


=  exp 


In 


JL  =  expflnlx  -  ju\ 
(TV  2  J 


A 

2 


ln(<j)- 


ln(2)^ 


where  max  =  highest  representable  number  in  binary  double  (64-bit)  float¬ 
ing-point  precision  (approximately  i.8xio3°8). 


The  absolute  value  is  computed  to  avoid  the  potential,  impossible  compu¬ 
tation  of  a  negative  logarithm.  The  negative  of  the  computed  term 


X- jU 


V2 


is 


inputted  into  the  error  function  for  x  <0<=>x-//<0. 

<tv2 


The  term 


jx-ju) 

a/2 


cannot  underflow  for  all  representable  x  and  fj.  since: 


— =-  —  iinii 

V2 

in  binary  double  (64-bit)  floating-point  precision. 


Computing  the  rearrangement  for  when  the  term  \x-ju\  overflows  allows 

to  compute  to  as  low  as  -L  rather  than 

a/2 


the  final  value  of  the  term 


X- JU 


2 


returning  00,  depending  on  the  value  of  the  term  a  in  the  denominator. 


A  numerically  stable  algorithm  of  the  error  function  for  any  real  number 
input  and  with  a  range  [-1,1]  computes  the  regularized  incomplete  gamma 

function  p(s,x)=  H^)=  '  r  ,.^e-  dr ;  using  either  the  series  or  continued 

r>)  r(s)J0 

fraction  representation: 


j-/>(i42)for  t<  0 
[ /'(■(,  ?2}  otherwise 
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Since  the  algorithm  for  computing  the  lower  incomplete  gamma  function 
is  only  stable  for  finite  representable  inputs,  where  t2  <  max ,  which  already 
computes  to  l  at  r  =  10,  the  error  function  is  set  to  -l  and  +1  for  t  <  — VTo  and 
t  >  vTo ,  respectively. 


Lognormal  distribution 

Conversion  of  mean  and  variance  to/from  ij  and  o 

Given  mean  and  variance,  it  is  possible  to  compute  the  unique  parameters 
for  the  lognormal  distribution  //  and  a  by  the  following  formulae: 


ju  =  In  (mean)  —  In  1 
2  v 


var 


mean" 


a  =  Jin  1  + 


var 

mean^ 


However,  both  of  these  expressions  are  susceptible  to  premature  overflow 
or  underflow  by  the  term  var  - .  The  formulas  may  be  made  numerically 

mean2 

stable  by  rearrangement. 


Firstly,  rearranging  the  term  —  -  to  the  term -  eliminates  the 

mean2  mean  V  mean  ) 

potential  of  premature  underflow  or  overflow  of  mean2,  thereby  modifying 
the  expression  for  q  and  <7  to  the  following: 


In  the  case  where  the  term -  overflows  the  expression  for  ji  and 

mean  V  mean  ) 

a  reduces  to  the  following  computable  forms  in  binary  double  (64 -bit) 
floating-point  precision: 

// =  ln(mean)-  — Inf  1h — Var  1 - mean  [mean] - ^  -  ln(mean)  - — Inf  Var  1  =  2  ln(mean)  -  —  ln(  var) 

2  v  mean  J  2  ^mean  J  2 


(7  =  Jin  1 


mean" 


-  1  [  var  | 

tor -  -  >max 

mean  V  mean) _ ^ ^  _  h^l 


mean 


=  ln(  var)  -  2  ln(mean) 
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Assuming  round-off  errors  are  negligible,  premature  overflow  is  eliminated  when 

calculating  //  by  computing  ju2  when  the  term  1 — [  var  j  overflows  and  compu- 

mean  1,  mean  J 

ting //i  otherwise: 

f  1  f  var  i 

_  I  n2 ,  tor -  -  >  max 

t*  mean  V  mean  ) 

/u{ ,  otherwise 

Similarly,  a  numerically  stable  algorithm  for  calculating  o  computes  a2 
when  the  term  — — f  -UU]  overflows  and  computes  oi  otherwise: 

mean  1  mean  ) 


a  = 


1  f  var 

<7  2 ,  ror -  -  >  max 

mean  (  mean  ) 

cr, ,  otherwise 


Underflow  of  the  term - 1  can  be  ignored  when  computing  fi  and  o 

mean  l  mean  ' 


since: 


1  +  0=1 


for  a  <  min  in  binary  double  (64-bit)  floating-point  precision,  where  min  = 
smallest  representable  number  in  binary  double  (64-bit)  floating-point 
precision  (approximately  4.9  x  10  324). 


Conversely,  it  is  possible  to  compute  mean  and  variance,  given  fj.  and  cr,  by 
the  following  formulae: 


cr 

JU  + - 

2 


mean  =  expj^ 
var  =  (expfcr2 )-  l)exp(2  /u  +  cr1 ) 


A  numerically  stable  algorithm  for  performing  these  operations  compute 
different  expressions  that  are  manipulated  yet  equivalent  for  different 
combinations  of  q  and  a. 


To  stably  compute  the  mean,  only  a  simple  change  in  the  order  of  opera¬ 
tions  is  necessary: 


mean  =  exp  fi 


cr 


~  exp 


/u  +  cr\ 


where  the  term,  oj  U  j  ?  overflows  only  when  the  true  value  of  the  mean 
overflows. 
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Two  rearranged  but  mathematically  equivalent  expressions  are  also  used 
in  different  situations  to  compute  variance: 


varj  =  (exp(cr)-l)  exp(2//  +  a1  )=  exp(ln[(exp(<r2 )-  l)ex  p(2//  + 
1 


(7~ 


=  exp 


In 


exp(<r2 ) 


+  2/j  +  2aa 


var,  =exp 


In 


exp  (<t  2 ) 


+  2  /i  +  2aa 


=  exp 


J 


In 


exi 


Y\ 


P(0"2  )J  2 

— - -  +  /U  +  CJ 


V  v 


Firstly,  the  variance  is  always  computed  to  zero  when  exp(<r2)is  equal  to  l. 
Var2  is  computed  whenever  2 [i  <  -max  and  2<7<7  >  max.  Vari  is  computed  in 
all  other  cases: 


var  =  < 


0,forexp(<r2)=  1 

var, ,  for  2/u  <  -  max  and  2  aa  >  max 
var, ,  otherwise 


The  variance  is  always  greater  than  the  largest,  positive  representable 
number  in  binary  double  (64-bit)  floating-point  precision  whenever  the 
term  <7 2  overflows. 


Whenever  either  2/u  or  200  overflows  but  not  both,  vari  =  var2  in  binary 
double  (64-bit)  floating-point  precision.  Underflow  of  the  term  in 

meani  and  the  term  2 00  in  vari  is  not  an  issue  since: 


scp(x  +  min)  =  exp(x) 


for  all  x  in  binary  double  (64-bit)  floating-point  precision.  The  term 

1 


In 


1- 


(exp(<r))‘T  _ 


could  never  underflow  since: 


[ln(x)|  >  min 

for  allx  7^  1  in  binary  double  (64 -bit)  floating-point  precision.  When  x  =  1, 
ln(x)  is  exactly  zero. 
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Computation  of  the  pdf 


As  with  the  Gaussian  distribution,  the  following  simple  formulation  of  the 
lognormal  pdf,  defined  for x>0,  ^  <=  9? ,  and  o  >  o,  is  also  susceptible  to 


premature  overflow  and  underflow  issues,  including  premature  overflow 
by  the  exponent  term  which  would,  in  turn,  cause  premature 

(ln(x)-//)2 


2o 


underflow  of  the  exponentiated  expression  exp 


2cr 


pdf(x) 


XG 


42n 


exp 


(ln(x)-  //) 

2<j2 


,2  "N 


A  simple  reformulation,  however,  eliminates  the  risk  of  premature  over¬ 
flow  and  underflow: 


pdf,  (x)  =  exp(ln(pdf  (x)))  = 


exp 


vlnl 


XG 


42jt 


exp 


(ln(x)-  //)' 


2 


2<t2 


z  exP 


^  +  21n(x)+21n(cr)+ln(2;r) 


cr 


where  the  term  ( ln^  F  )  overflows  only  when  the  true  value  of  the  pdf  is 


zero. 


Underflow  of  the  term  f  ln(x)-/d|  js  not  an  issue  since 


<7 


exp 


-—(a  +  min) 
2  v 


=  1  for  a  e  {F :  a  ^  a  ±  min} 


where  F  is  the  set  of  all  numbers  representable  in  binary  double  (64-bit) 
floating-point  precision  and  a  is  set  as  the  expression,  2in(x)+2in(o-)+in(2^). 


Since  the  above  rearrangement  always  computes  NaN  for  x  =  o,  the  pdf  of 
x  =  min  is  computed  instead  in  this  case,  which  may  or  may  not  necessari¬ 
ly  limit  to  zero  depending  on  the  values  of  p  and  o.  Thus,  the  final  algo¬ 
rithm  for  computing  the  pdf  is  as  follows: 


pdf(x)  = 


(pdf  (min),  for  x  =  0 
[  pdf,  (x),  otherwise 
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Computation  of  the  cdf 

Where  the  error  function  is  numerically  stable  for  any  real  number  input 
and  has  the  range  [-1,1],  the  following  formula  for  the  lognormal  cdf  is  ac¬ 
curate  only  when  computation  of  the  term  is  numerically  stable: 

crV2 

A  numerically  stable  calculation  of  the  term  is  done  simply  by 

crV2 

computing  parts  of  it  in  a  certain  order  as  follows: 

ln(x)-//  _  I  ( ln(x)-/A 

<rV  2  a  {  J 

Here,  the  term  [  ln^b  £  j  can  never  overflow  or  underflow  by  itself  for  all 

finite  p  and  nonzero  x  in  binary  double  (64-bit)  floating-point  precision 
since: 


min 

— =  min 

V2 


and 


ln(x)  +  max  =  +  max 

for  all  finite  x  in  binary  double  (64 -bit)  floating-point  precision. 

The  error  function  is  computed  with  the  same  algorithm  described  for 
computation  of  the  Gaussian  cdf. 

Exponential  distribution 

Parametrization 

The  exponential  distribution  is  characterized  by  a  single,  scale  parameter  (3 
that  is  equivalent  to  the  mean.  The  variance  is  the  square  of  the  mean. 
Since  conversions  are  one  step,  no  rearrangement  is  necessary  or  possible: 


mean  =  J3 
var  =  / 3 2 
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Computation  of  the  pdf 

Without  computing  a  rearranged  form,  the  following  expression  for  the 
exponential  pdf,  defined  for  x  >  o  and  /3  >  O,  could  result  in  premature  un¬ 
derflow  of  the  exponentiated  expression,  exp|--hj,  when  /3  <  l: 

pdf  (x)  =  —  expf-  — 

/?  \  P) 

The  following  algebraically  equivalent  but  rearranged  expression  is  accu¬ 
rate  for  all  x  and  / 3 : 


pdf(x)  =  -^-exp] 


'  x^ 

=  exp 

in 

f 1  1 

-rexp 

'  x^ 

| 

=  exp 

-^-ln  (P) 

,  Pj 

l  P  ' 

v  P) 

1  P  \ 

where  the  term  T  overflows  only  when  the  true  value  of  the  pdf  is  zero. 


Underflow  of  the  term  —  is  not  an  issue  since: 

P 


exp[a  -  min]  =  1  for  aejf  :«/«  +  mir 


where  F  is  the  set  of  all  numbers  representable  in  binary  double  (64-bit) 
floating-point  precision  and  a  is  set  as  the  expression  -ln(/?). 


Computation  of  the  cdf 

The  exponential  cdf  can  be  stably  computed  by  the  following,  simple 
expression  for  all  x  and  /3 : 


cdf(x)  =  1  -  exp|  -  — 


Computation  of  the  quantile  (cdf1) 

Since  the  cdf  expression  can  be  easily  solved  for  x,  a  closed-form  formula 
for  the  exponential  quantile  exists,  which  is  stable  for  all  x  and  /?: 

x  =  -J3  ln(l  -  cdf  (x)) 
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Gamma  distribution 

Conversion  of  mean  and  variance  to/from  a  and  /3 

Given  mean  and  variance,  it  is  possible  to  compute  the  unique  parameters 
for  the  gamma  distribution,  a  and  jS,  by  the  following  formulae: 


mean' 

a  = - 

var 


P 


var 

mean 


While  the  expression  for  (3  is  stable  for  all  mean  and  variance,  a  should  be 
computed  as  follows: 


mean2  .  mean 

- ,  tor - >  max 

var  var 

mean^l  , 

-  mean,  otherwise 

var  J 


The  converse  operation  computes  mean  and  variances  of  the  gamma  dis¬ 
tribution  from  a  and  (5,  by  the  following  formulae: 


mean  =  a[3 
var  =  a/3 2 

The  expression  for  mean  is  stable  for  all  a  and  (3,  but  variance  should  be 
computed  as  follows: 


var  =  a/32  =  /3{a/3 ) 


Computation  of  pdf 

The  following  expression  for  the  gamma  pdf,  defined  for  x  >  o,  a  >  O,  and 
P  >  O,  must  be  formulated  carefully  to  prevent  inaccuracies  from  numer¬ 
ous  potential  premature  overflow  and  underflow  issues: 


pdf(x)  = 


r  (a)/r 


-exp 


But,  to  obtain  a  stable  rearrangement  of  the  gamma  pdf,  it  is  first  neces¬ 
sary  to  notice  how  the  term  r(«)  is  computed  using  the  following  expres¬ 
sion  for  Lanczos  approximation  of  the  gamma  function  (Lanczos  1964): 
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(  1  y+2 

f  n 

■  a  +  y  +  —  exp 

-  \a+r+-\ 

V  2) 

L  v  2  jj 

aM) 


with 


Ar  («) = ir)  +  a  ir  +  Pi  ir ) 

2  a  + 1 


a 


(a-l) 


(a  +  l\a  +  2) 


and 


pAr)=Xc% 

<5=0 


+  y  + 


+  /  + 


where  y  is  an  arbitrarily  chosen  constant  such  that  +  y  +  ^ j  >  0  so  that 

the  infinite  series,  Ar(a),  converges  as  required  and  c)  denotes  the  coeffi¬ 
cient  of  the  j-th-degree  term  in  the  i-th-degree  Chebyshev  polynomial  of 
the  first  kind. 


Since  r(a  +  i)  =  ar(«),  we  have: 


r(a)=lk±l) 

a 


i 


r  1 

(  n 

a  +  y +  —  exp 

-\a  +  y  +  -\ 

1  2) 

L  V  2j\ 

Since  the  series,  Ar(a),  is  convergent,  it  may  be  truncated  to  obtain  an  ap¬ 
proximation  with  desired  precision.  For  fixed  y  and  truncation  order  N, 
computational  efficiency  is  enhanced  by  resolving  the  rational  fractions  in 
the  series  Ar(a)  into  their  constituent  partial  fractions  so  that: 


Ma) 


a 


yA,(/)  +  Pi(/)  , 

2  a+ 1 


+  Pi(y) 


t(a- 1) 


(i a  +  \\a  +  2 ) 


'K  (r)+ 


y  b,Xr) 

a  +  n 


where  the  coefficients  b(y)  are  pre-computed  for  fixed  y,  which  are  related 
to  the  original  coefficients  p(r)  of  the  series  Ar(a)  and  similarly  indepen¬ 
dent  of  a. 


To  illustrate  how  the  coefficients  b(y)  are  pre-computed  after  decomposing 
rational  fractions  of  the  series  Ar(a)  the  process  is  shown  for  a  simpler 

case,  when  the  truncation  order  is  3.  First  the  following  is  obtained  by 
combining  like  terms  after  partial-fraction  decomposition  of  the  rational 
fractions  in  the  truncated  series: 
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Ar  (r)  ~  v  a  (r ) + a  (r)— - + a  (r)j  ^  ,  +  a  ir),  a}°L  ^  , 

\a  +  \\a  +  2)  [a  + 1  \a  +  2  \a  +  3) 


a  + 1 


-  b.M+ aOV aM+  p,(rh  zPkhMrtlEM + -_*pM±^pM  +  -J^pM 


a  +  2  a  +  3 


a +  2 


a  +  3 


Then  we  notice  that: 


4  to  -  i  A  (rH AH -aOAAAA) + zAfeHjfA) + AA> 

2  a  +  1  a  +  2  a  +  3 

=„„W+Mzl+Mrl+Mi=(>„w4'>.W 


a  +  1  a  +  2  a  +  3 


“  a  +  « 


when  we  let: 


(r) = |  a  (r) + a  (r) + a  (r) + a  (r) 
*i(r)= -a  (r) + 2  a  (r)  -  3  a  M 

^2(r)=-6p2(r)+24p3(r) 

*3(/)  =  -30a(/) 


After  examining  the  dependence  of  relative  error  as  a  function  of  y  and  the 
truncation  order  of  the  series,  Ar(a),  it  has  been  determined  that  using  y  = 
10.900511  and  a  truncation  order  of  10  guarantees  16-digit  floating-point 
accuracy  required  for  implementation  in  binary  double  (64-bit)  floating¬ 
point  precision,  which  yields  the  following  coefficients  (Pugh  2004): 


b0  =  2.48574089138753565546 xlO-5 
b ,  =1.05142378581721974210 
b2  =-3.45687097222016235469 
b3  =4.51227709466894823700 
b4  =-2.98285225323576655721 
b5  =  1.05639711577126713077 
b6  =  -1.95428773191645869583  x  10”1 
b7  =  1.70970543404441224307 xl0“2 
bs  =  -5.719261 17404305781283  x  10“4 
b9  =  4.63399473359905636708 xlO-6 
bl0  =  -2.71994908488607703910x  10-9 


Currently,  EASEE  uses  an  implementation  of  the  Lanczos  approximation 
by  M.  T.  Flanagan  in  his  Java  Scientific  and  Numerical  Library  with  y  = 
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5.0,  a  truncation  order  of  6,  and  an  error  bound  less  than  2xio10,  which 
has  the  following  coefficients: 

b0  =1.000000000190015 
bt  =76.18009172947146 
b2  =-86.50532032941677 
b3  =24.01409824083091 
b4  =-1.231739572450155 
b5  =  0.1208650973866179xl0'2 
b6  =  -0.5395239384953x  10'5 

Now  return  to  the  original,  problematic  formulation  of  the  gamma  pdf: 


pdf(x)  = 


r  (a)p‘ 


-exp 


x 

~P 


and  set: 


r(a)  = 


■sflft 


(  1 Y+2 

(  1Y 

-\cc  +  y  +  -\  exp 

_T+r+2j_ 

(  1  x+2 

f  1) 

-  a  +  y  +  —  exp 

-  «+r+- 

V  2  J 

L  v  2J] 

A(a) 


bM +i^ 

a +  n 


to  obtain: 


pdf  (x)  = 


Y{a)p‘ 


-exp 


'  x ' 


P 


x“  1  expi 


^  x^ 


v  Pj 


V  H  J 


■Jin  ( 


a  \ 


a  +  y  +  - 


exp 


-]  a  +  r+- 


b0(r)  +  Yj 


b„(r ) 


,  a +  n 


P° 


which  is  the  actual,  full  expression  of  the  approximated  gamma  pdf. 

This  formulation  of  pdf  is  unstable  when  computed  within  binary  double 
(64-bit)  floating-point  precision,  mostly  due  to  premature  overflow  by  the 
approximation  expression  for  the  term  r(«)  and  the  term  x“"‘  as  well  as 

premature  underflow  by  the  term  exp|-  ~)  .  In  fact,  there  are  many  (even 

relatively  small)  values  for  x  and  a  that  cause  r(«)  and  x“ 1  to  overflow  si¬ 
multaneously,  since  both  are  rapidly  increasing  functions  as  a  increases, 
generating  NaN  due  to  the  computation  of  00/00,  even  though  the  true  ratio 
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of  - —  is  a  simple,  finite  number  that  is  easily  representable  in  binary 

r(«) 

double  (64-bit)  floating-point  precision.  NaN  may  also  be  returned  even 
when  only  x“_1  overflows,  if  the  term  6XP(“)  underflows  to  zero,  which 


results  in  the  indeterminate  form,  Oxoo,  in  the  numerator.  These  observa¬ 
tions  along  with  many  other  potential  scenarios  that  could  cause  inaccu¬ 
rate  (or  more  often  unviable)  solutions  to  the  gamma  pdf  make  this  formu¬ 
lation  quite  unpredictable  and  unsatisfactory. 


Firstly,  a  few  basic  manipulations  of  the  original  formula  result  in  the 
following,  more  stable  form: 


Pdf,(x)  = 


T(a)p‘ 


-exp 


xa  1  exp] 


p2jr  ( 
a 


a  +  y  +  - 


1 

nC(  H — 
2 


exp 


—  a  +  y  +  - 


“  a  +  n 


Pa 


axa  1  exp]  a  +  y  + 


(  „=1  a  +  n)\ 


exp 


:  exp 


:  exp 


In 


axa  1  exp]  a  +  y  + 


/rV2^k(r)+Z—  II  «  +  r  + 

(  „=i  a +  n  " 


A-  —  +  aB 

P 


1  xa+-  / 

n  2  1 


exp 


x 

kPjj 


In  a  +  y  + 


and 


where  ^4  =  ln(a)-ln(x)+a  +  ^  +  —  (y)  +  V  ^ 

2  2  (  a  +  n  j 

b  =  in(x)-  in(/?)-  ln^a  +  y  +  ,  which  is  stable  for  most  x,  a,  and  /?.  When  the 

term  A  overflows  and  B  >  o,  the  following  is  computed  instead: 


pdf2(x)  = 


r  («)/?“ 


(  X  ^ 

-expj 

V~pj 

| «  exp 

A- 

— h  aB 

p 

=  exp 

(Ap~x  +  ( a/3)B ) 


P 
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To  avoid  automatically  computing  a  pdf  value  of  +oo  for  x  =  o  using  these 
rearrangements,  the  pdf  of  x  =  min  is  returned  instead  as  a  very  close  ap¬ 
proximation,  which  would  only  limit  to  zero  for  sufficiently  large  a. 

Thus,  the  final  algorithm  for  the  gamma  pdf  computation  is  as  follows: 

pdf  (min),  forx  =  0 
pdf,  (x),  for  T  >  max  and  B  >  0 
pdf[  (x)  otherwise 

Here,  the  term  A  never  overflows  for  all  nonzero,  finite,  positive  values  of  x 
and  a  in  binary  double  (64-bit)  floating-point  precision  and  y  =  5.0,  as  re¬ 
quired  to  prevent  premature  overflow. 

To  confirm  this,  we  first  notice  that  the  expressions  60(/)+ and 

n=\  Ot  ~f~  Tl 

«  +  r+y  are  nonzero  and  finite  in  binary  double  (64-bit)  floating-point  pre¬ 
cision  in  the  given  domain  of  x,  a,  and  y  (i.e.,  x  >  o,  a  >  O,  and  y  =  5.0). 
This  is  easy  to  see  with  the  expression  a  +  r  +  y  =  «  +  5  +  y  which  can  never 

equal  zero  because  it  is  the  summation  of  nonzero,  positive  numbers  and 
can  never  overflow  since: 

max+  5  +  —  =  max 
2 

in  binary  double  (64-bit)  floating-point  precision.  The  expression 

N  b  (y)  N  b  (5) 

^o(z)+y^iL^  =  ^o(5)+y^ili-  can  never  overflow  by  examination  ofthe  bis) 

~la  +  n  a +  n 

coefficients  for  truncation  order  6  and  noting  that  a  +  n  >  1  for  any  positive, 
nonzero  integer,  n.  It  can  never  equal  zero  since  b0(y)+Y  =  0  implies 

n=\  OC  +  Tl 

that  r(«)  =  0 ,  which  contradicts  the  fact  that  r(«)  >  0  for  a  >  0 . 

Hence,  given  that: 

-744.4400719213812  ~  ln(min)<  ln(a)<  ln(max)«  709.782712893384  for  ae{F:a*  0}- 

where  F  is  the  set  of  all  numbers  representable  in  binary  double  (64-bit) 
floating-point  precision,  we  have: 
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,  .  .  ,  v,  lna  +  )'  +  - 

A  -  In  (a)  - 1  n  (x )  +  a  +  y  +  —  -  ~K  -  ln(  b0  (y)  +  V  ^  | - —  <  a  +  f  +  —  +  5  x  |ln(min)| 1 

22  (  ~l  a  +  nj  2  2 


The  expression  a  +  r  +  ^  +  5x |in(min)|  =  «  +  5  +  ^  +  5 x]in(min)|  can  never  overflow 


since: 


max+  5  +  i  +  5  x  |ln(min)|  =  max 

in  binary  double  (64-bit)  floating-point  precision. 


The  term,  B,  never  overflows  for  all  nonzero,  finite,  positive  values  of  x,  a, 
and  /?  in  binary  double  (64-bit)  floating-point  precision  and  y  =  5.0: 


\B  = 


ln(jc)  —  ln(/?)-  Inj"  a  +  y  +  ^  j 


<  3|ln(min)| 


<  max 


Underflow  of  ^4  is  not  an  issue  since: 


max 

in  binary  double  (64-bit)  floating-point  precision. 

Underflow  of  any  of  one  or  more  of  the  terms  A,  U,  and  aB  may  be  ignored 
since: 

exp(a  ±  min)  =  1  for  a  e  {F  :  a  a  ±  min} , 
exp(a  +  2 min)  =  1  for  a  e  {F  :  a  *  a  +  2 min},  and 
exp(+  3  min)  =  1 

where  F  is  the  set  of  all  numbers  representable  in  binary  double  (64-bit) 
floating-point  precision  and  a  is  set  as  the  expression  for  the  non¬ 
underflowing  term(s). 

Computation  of  the  cdf 

Like  the  error  function  in  the  Gaussian  and  lognormal  cdfs,  the  gamma  cdf 
is  also  related  to  the  regularized  incomplete  gamma  function 
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P{s,x)  = 


r{s,x ) 

r(0 


rs~'e~'  dr ,  which  is  likewise  computed  using  either  the  se¬ 


ries  or  continued  fraction  representation: 


cdf(x)  = 


f  x 

a , — 

V  P 


r(«) 


While  the  individual  terms  yja,Tj  and  r(a)  overflow  at  modest  values  of  x, 

a,  and  (3,  the  regularized  incomplete  gamma  function  itself  never  overflows 
since  it  has  the  limiting  values: 


p(s,o)  =  o  and  p(s,°o)=i. 


Thus,  it  is  better  to  implement  the  exponentiation  of  the  difference  of  the 
logarithms  of  the  two  terms: 


(  \  rr’7 

cdfW=p(«^]=A^ 


=  exp 


In 


r  a 


p 

r(a) 


:  exp 


(  ( 


In 


a, — 

v  P 


-  ln(r(a)) 


With  this  rearrangement,  the  gamma  cdf  is  computable  for  a  much  more 
expansive  domain  of  x,  a,  and  /?,  since  it  is  only  required  that  the  loga¬ 
rithm  of  the  two  terms  j  and  r(«)  not  overflow.  For  the  vast  majority 

of  cases,  this  rearranged  formulation  of  the  gamma  cdf  is  stable.  In  the 
rare  instance  that  the  rearrangement  proves  unstable,  the  algorithm  com¬ 
putes  a  direct  numerical  integration  of  the  gamma  pdf,  which  is  always 
stable.  Information  on  the  implementation  of  the  Gauss-Legendre  quadra¬ 
ture  rule  for  efficiently  integrating  pdfs  is  given  in  the  final  section  of  the 
report  on  numerical  issues. 


X  — *  transformed  Rice-Nakagami  distribution 

The  distributions  of  the  powers  of  certain  acoustic  and  seismic  signals  are 
closely  approximated  by  the  Rice-Nakagami  model.  Specifically,  the  sig¬ 
nal-power  distribution  resembles  the  X  X1  transformed  Rice-Nakagami 
distribution  when  the  amplitude  is  Rice-Nakagami  distributed,  since  sig¬ 
nal  power  is  the  square  of  signal  amplitude.  For  consistency,  all  sensors 
are  modeled  to  work  with  signal  power  rather  than  amplitude,  thus  the 
X  X2  transformed  Rice-Nakagami  distribution  is  encoded.  It  turns  out 
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that  it  is  also  simpler  to  implement  the  X  -» X1  transformed  Rice- 
Nakagami  distribution,  whose  expressions  for  computing  mean  and  va¬ 
riance  from  its  unique  parameters  v  and  a  are  clean  polynomials  and,  un¬ 
like  the  untransformed  version,  do  not  involve  the  half-degree  Laguerre 
polynomial. 

Derivation  of  the  X  — ►  X2  transformed  Rice-Nakagami  distribution 

Statistically,  the  x  ->  x2  transformed  Rice-Nakagami  distribution  is  the 
distribution  of  X2  where  X  is  a  random  variable  with  a  Rice-Nakagami  dis¬ 
tribution.  A  fundamental  theorem  for  random  variable  transformation 
states  that,  if  pT(t)  is  the  value  of  the  probability  density  of  the  continuous 
random  variable  T  at  t  and  the  function  x  =  ® (t)  is  differentiable  and  mono¬ 
tonic  for  all  values  within  the  range  of  T  for  which  its  probability  density 
does  not  equal  o  (i.e.,  pT{t)  *  o),  the  probability  density  of  X  is  given  by: 

r.  ./  VI  d&'Hx)  . -d0(t )  n 

pMApM 

0,  otherwise 

Hence,  for  the  Rice-Nakagami  distribution,  where  x  =  ®(?)  =  t2 ,  we  have: 

0,  otherwise 

where: 

<t>~'{x)  =  +4x  since  t>o, 
d®-'{x)  =  d[jx)=  i  .  and 

dx  dx  2  4x 

dt 
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Thus,  for  t  *  o: 


p  X{x)=  Pt[4x- 


4% 

— fexp| 

<T 

1 

~~  2cr2 


-jx+v2) 

2a2 

-{x  +  v 
2<r2 


2y[x\ 


v-2  V 


1 

2-v/x 


Vyfx 


where  I0  is  the  zeroth  order  modified  Bessel  function  of  the  first  kind  and  v 
and  <7  are  the  nonnegative  parameters  of  the  Rice-Nakagami  distribution. 


Parameters  of  the  X—*  X2  transformed  Rice-Nakagami  distribution 

By  definition,  the  mean  of  transformed  pdf,  px(x),  is  the  integral  of  all  val¬ 
ues  within  the  range  of  X  with  respect  to  (i.e.,  weighted  by)  its  probability 
density: 


00 

E(x)=  [  xpx(x)dx 

—oo 


which,  for  .r  =  ® {t)=t2  and  t  +  o,  is  equal  to  the  second  raw  moment  p2  of 
the  untransformed  pdf  pT{t): 


E(x)=  j '  xpx(x)dx  =  | <i>(t)pr [<Z>  ‘(x)  • 


d>P^{x) 


dx 


o° 

dx=  j  (t)dt  =  E(o(t))  =  e(t2)=  p2 


! 

Since  //2  =  2<j2  +  v2  for  the  Rice-Nakagami  distribution,  we  have: 

E{x)  =  p’_  =2a2  +v2 

Also  by  definition,  the  variance  of  the  transformed  random  variable,  X,  can 
be  expanded  as  follows: 

Var{x)  =  e[{X  -  E{x))2  ]  =  e[x2  -  2 XE{x)  +  E{x)2  ] 

=  E(X2)-2E{X)E{X)  +  E{X)2  =e(x2)-2E{x)2  +  E{x)2  =  e(x2)-E{x)2 

which,  for  t  =  ®(r)  =  r  and  t  +  o,  may  be  expressed  using  the  second  and 
fourth  raw  moments  and  of  the  untransformed  pdf  pT(t): 


Var{x)=E[x2)-E{X )2  =  e(t4)- e(t2J  =  -^2'j  =  (8<r4  +  8ctV  +  v4)-  {2a2  +  v2f 

=  (8<t4  +  8 a2v2  +  v4)-  (4<r4  +  4cr2v2  +  v4)=  8<r4  +  8<t2|'2  +  v4  -  4<r4  -  4<j2v2  -  v4  =  4<r4  +  4<t21'2 
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Given  e(x)  =  2<j2  +  wand  Var(x)=4a4  +4<rV  for  the  .v ->x2  transformed  Rice- 
Nakagami  distribution,  the  parameters,  a  and  v,  of  the  transformed  pdf 
px(x)  maybe  obtained  from  £’(x)and  Var(x  )by  the  following  formulae: 


a  = 


E(x) -  sj E(x )2  -  Var(x) 


v  =  [e(x)2  -  Fflr(z)}1 


Conversion  of  mean  and  variance  to/from  v  and  o 


From  the  previous  section,  we  have  the  following  formulae  for  computing 
the  parameters,  <7  and  v,  from  mean  and  variance  for  the  X  X2  trans¬ 
formed  Rice-Nakagami  distribution: 


mean  -  Vmean2  -  var 

V  2 

v  =  [mean2  -  var]-* 

Conversely,  we  have  the  following  expressions  for  mean  and  variance  with 
respect  to  the  parameters,  a  and  v,  for  the  x^x2  transformed  Rice- 
Nakagami  distribution: 


mean  =  2cr2  +v 2 
var  =  4<r4  +  4a2v2 


While  the  exact  details  are  not  included  here,  these  conversion  expressions 
are  made  stable  using  similar  techniques  used  for  the  other  previously  de¬ 
scribed  distributions. 

Computation  of  the  pdf 

As  derived  earlier,  the  X  ->  X2  transformed  Rice-Nakagami  pdf,  defined  for 
x >  o ,  v>o,  and  o- > o ,  is  as  follows: 

rwTi 

where  I0  is  the  zeroth  order  modified  Bessel  function  of  the  first  kind. 
Again,  details  on  implementing  this  pdf  are  not  included  here.  As  with  the 
gamma  function  in  the  gamma  pdf,  it  is  important  to  examine  and  dissect 
the  approximating  method  used  to  compute  I0  when  manipulating  the  en¬ 
tire  pdf  expression  into  stable  rearrangements. 


pdf(x)  = 


2(7 


-exp 


2rr2 
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Computation  of  the  cdf 

To  compute  the  X  ->  X2  transformed  Rice-Nakagami  cdf,  the  transformed 
pdf  is  numerically  integrated  using  Gauss-Legendre  quadrature.  While  this 
approach  is  generally  more  computationally  intensive  than  a  direct  ap¬ 
proximation,  its  stability  is  guaranteed  if  the  algorithm  for  the  pdf  compu¬ 
tation  is  fully  stable.  More  importantly,  it  is  a  simpler  alternative  to  devis¬ 
ing  a  stable  approximation  of  the  first-order  Marcum  Q-Function  in  the 
X  -» X2  transformed  Rice-Nakagami  cdf.  An  overview  on  implementing 
Gauss-Legendre  quadrature  to  numerically  integrate  pdfs  is  in  the  final 
section  of  this  appendix. 

As  with  the  pdf,  it  is  also  possible  to  derive  the  x->x2  transformed  Rice- 
Nakagami  cdf.  In  general,  any  continuous  random  variable  X  that  is  re¬ 
lated  to  another  continuous  random  variable  T  by  a  function  x  =  ®(r)  has 
the  following  cdf: 


Fx {x)=P{X<x)  =  P{0{T)<x)=  P{{t  e  V: 0{t)<x})=  J  fT (t)dt 

where  W  is  the  support  set  of  T,  Fx  is  the  cdf  of  X,  and/r  is  the  pdf  of  T. 

Given  that  the  function  0  is  strictly  increasing,  as  is  the  case  for 
x  =  ®(r)  =  t2,  we  have: 


0  '(x) 

Fx(x)=  J  fT(t)dt  =  p{t<0-1(x))=  \fT(t)dt  =  FT{0-'(x)) 

{t&:0(T)<x}  -to 


where  Ft  is  the  cdf  of  T. 


Since  the  untransformed  Rice-Nakagami  distribution  has  the  following 
cdf: 


cdf(x)=l 

we  have  the  following  expression  for  the  X  -» X2  transformed  Rice- 
Nakagami  cdf: 


Fx(x)=G(^-‘(x))=i-a 


0-'{x) 


=i -a 


V  a/x 

< j  ’  a 
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where  <? {x)=+4x ,  since  =  {t  e  9? :  t  >  0}  and  Qi  is  the  first-order  Marcum  Q- 
Function.  A  numerical  integration  of  the  pdf  is  favored  in  this  case  to  avoid 
formulating  stable  implementations  of  Qi,  the  first-order  Marcum  Q- 
Function. 


Stable  implementation  of  Gauss-Legendre  quadrature  on  arbitrary  integration 
intervals 

Given  the  abscissas  {xNk  and  weights  {wN  k  for  the  AT-point  Gauss- 

Legendre  rule  over  [-1,1],  it  is  possible  to  apply  the  rule  on  a  function  f(t ) 
over  an  arbitrary  interval  [a, b]  by  the  following  change  of  variable: 

t  =  a1±+b-ax  and  dtJ_ZdLdx 
2  2  2 


Then  the  following  relationship: 

is  used  to  obtain  this  AT-point  Gauss-Legendre  quadrature  formula: 


b  -  a 
2 


b  —  a 
2 


XN,k 


For  the  technique  to  stably  integrate  a  function  f(t)  over  the  interval  [ a,b ], 
it  is  only  required  that  the  abovementioned  change-of-variable  operation 
is  made  stable,  provided  that  the  algorithm  for  the  AT-point  Gauss- 
Legendre  rule  stably  computes  all  the  abscissas  and  weights  over  [-1,1]  for 
some  specified  N  and  that  the  function  f(t )  is  stable  over  the  specified 

interval  [ a,b ].  For  integrating  pdfs  with  nonnegative  support,  the  expres¬ 
sion  (b  -  a)  /  2  is  always  stable.  However,  whenever  the  expression 
(a  +  6)  /  2  overflows,  it  should  be  computed  using  the  distributed  formula 
x  /  2  +  y  /  2  (but  only  then).  Thus,  we  have  for  all  pdf  with  nonnegative 
support: 


cdf  ( b )  =  |  pdf  (t)dt  = 


(b  -  a) 


Z,,  U  U 

wN,kf  y  +  y  + 


( b  -  a) 


b  (b- a) 


2X,/ 


(a  +  b)  ( b-a ) 


|,  for  a  +  b>  max 
L  otherwise 


(a  +  b) 

Since  - - -  >  max  <=>  a  +  b>  max  • 
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